smm_objective

 Evaluates the objective of the Simulated Moments Method.

 INPUTS:
  xparams          [double]  p*1 vector of estimated parameters. 
  sample_moments   [double]  n*1 vector of sample moments (n>=p).
  weighting_matrix [double]  n*n symetric, positive definite matrix.
  options          [      ]  Structure defining options for SMM.
  parallel         [      ]  Structure defining the parallel mode settings (optional).

 OUTPUTS: 
  r                [double]  scalar, the value of the objective function.
  junk             [      ]  empty matrix.

 SPECIAL REQUIREMENTS
  The user has to provide a file where the moment conditions are defined.

0001 function [r,flag] = smm_objective(xparams,sample_moments,weighting_matrix,options,parallel)
0002 % Evaluates the objective of the Simulated Moments Method.
0003 %
0004 % INPUTS:
0005 %  xparams          [double]  p*1 vector of estimated parameters.
0006 %  sample_moments   [double]  n*1 vector of sample moments (n>=p).
0007 %  weighting_matrix [double]  n*n symetric, positive definite matrix.
0008 %  options          [      ]  Structure defining options for SMM.
0009 %  parallel         [      ]  Structure defining the parallel mode settings (optional).
0010 %
0011 % OUTPUTS:
0012 %  r                [double]  scalar, the value of the objective function.
0013 %  junk             [      ]  empty matrix.
0014 %
0015 % SPECIAL REQUIREMENTS
0016 %  The user has to provide a file where the moment conditions are defined.
0017 
0018 % Copyright (C) 2010-2012 Dynare Team
0019 %
0020 % This file is part of Dynare.
0021 %
0022 % Dynare is free software: you can redistribute it and/or modify
0023 % it under the terms of the GNU General Public License as published by
0024 % the Free Software Foundation, either version 3 of the License, or
0025 % (at your option) any later version.
0026 %
0027 % Dynare is distributed in the hope that it will be useful,
0028 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0029 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0030 % GNU General Public License for more details.
0031 %
0032 % You should have received a copy of the GNU General Public License
0033 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0034 
0035 global M_ options_ oo_
0036 persistent mainStream mainState
0037 persistent priorObjectiveValue
0038 
0039 flag = 1;
0040 
0041 if nargin<5
0042     if isempty(mainStream)
0043         if matlab_ver_less_than('7.12')
0044             mainStream = RandStream.getDefaultStream;
0045         else
0046             mainStream = RandStream.getGlobalStream;
0047         end
0048         mainState  = mainStream.State;
0049     else
0050         mainStream.State = mainState;
0051     end
0052 end
0053 
0054 if isempty(priorObjectiveValue)
0055     priorObjectiveValue = Inf;
0056 end
0057 
0058 penalty = 0;
0059 for i=1:options.estimated_parameters.nb
0060     if ~isnan(options.estimated_parameters.upper_bound(i)) && xparams(i)>options.estimated_parameters.upper_bound(i)
0061         penalty = penalty + (xparams(i)-options.estimated_parameters.upper_bound(i))^2;
0062     end
0063     if ~isnan(options.estimated_parameters.lower_bound(i)) && xparams(i)<options.estimated_parameters.lower_bound(i)
0064         penalty = penalty + (xparams(i)-options.estimated_parameters.lower_bound(i))^2;
0065     end
0066 end
0067 
0068 if penalty>0
0069     flag = 0;
0070     r = priorObjectiveValue + penalty; 
0071     return;
0072 end
0073 
0074 save('estimated_parameters.mat','xparams');
0075 
0076 % Check for local determinacy of the deterministic steady state.
0077 noprint = options_.noprint; options_.noprint = 1;
0078 [local_determinacy_and_stability,info] = check(M_,options_,oo_); options_.noprint = noprint;
0079 if ~local_determinacy_and_stability
0080     r = priorObjectiveValue * (1+info(2));
0081     flag = 0;
0082     return
0083 end
0084 
0085 simulated_moments = zeros(size(sample_moments));
0086 
0087 % Just to be sure that things don't mess up with persistent variables...
0088 clear perfect_foresight_simulation;
0089 
0090 if nargin<5
0091     for s = 1:options.number_of_simulated_sample
0092         time_series = extended_path([],options.simulated_sample_size,1);
0093         data = time_series(options.observed_variables_idx,options.burn_in_periods+1:options.simulated_sample_size);
0094         eval(['tmp = ' options.moments_file_name '(data);'])
0095         simulated_moments = simulated_moments + tmp;
0096         simulated_moments = simulated_moments / options.number_of_simulated_sample;
0097     end
0098 else% parallel mode.
0099     if ~isunix
0100         error('The parallel version of SMM estimation is not implemented for non unix platforms!')
0101     end
0102     job_number = 1;% Remark. First job is for the master.
0103     [Junk,hostname] = unix('hostname --fqdn');
0104     hostname = deblank(hostname);
0105     for i=1:length(parallel)
0106         machine = deblank(parallel(i).machine);
0107         if ~strcmpi(hostname,machine)
0108             % For the slaves on a remote computer.
0109             unix(['scp estimated_parameters.mat ' , parallel(i).login , '@' , machine , ':' parallel(i).folder ' > /dev/null']);
0110         else
0111             if ~strcmpi(pwd,parallel(i).folder)
0112                 % For the slaves on this computer but not in the same directory as the master.
0113                 unix(['cp estimated_parameters.mat ' , parallel(i).folder]);
0114             end
0115         end
0116         for j=1:parallel(i).number_of_jobs
0117             if (strcmpi(hostname,machine) && j>1) || ~strcmpi(hostname,machine)  
0118                 job_number = job_number + 1;
0119                 unix(['ssh -A ' parallel(i).login '@' machine ' ./call_matlab_session.sh job' int2str(job_number) '.m &']);
0120             end
0121         end
0122     end
0123     % Finally the Master do its job
0124     tStartMasterJob = clock;
0125     eval('job1;')
0126     tElapsedMasterJob = etime(clock, tStartMasterJob);
0127     TimeLimit = tElapsedMasterJob*1.2;
0128     % Master waits for the  slaves' output...
0129     tStart = clock;
0130     tElapsed = 0;
0131     while tElapsed<TimeLimit
0132         if ( length(dir('./intermediary_results_from_master_and_slaves/simulated_moments_slave_*.dat'))==job_number )
0133             break
0134         end
0135         tElapsed = etime(clock, tStart);
0136     end
0137     try
0138         tmp = zeros(length(sample_moments),1);
0139         for i=1:job_number
0140             simulated_moments = load(['./intermediary_results_from_master_and_slaves/simulated_moments_slave_' int2str(i) '.dat'],'-ascii');
0141             tmp = tmp + simulated_moments;
0142         end
0143         simulated_moments = tmp / job_number;        
0144     catch
0145         r = priorObjectiveValue*1.1;
0146         flag = 0;
0147         return
0148     end
0149 end
0150 
0151 r = transpose(simulated_moments-sample_moments)*weighting_matrix*(simulated_moments-sample_moments);
0152 priorObjectiveValue = r;
0153 
0154 if (options.optimization_routine>0) && exist('optimization_path.mat')
0155     load('optimization_path.mat');
0156     new_state = [ r; xparams];
0157     estimated_parameters_optimization_path = [ estimated_parameters_optimization_path , new_state ];
0158     save('optimization_path.mat','estimated_parameters_optimization_path');
0159 end