moment_function

 Evaluates the moment function of the Simulated Moments Method (discrepancy between sample and
 ).

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

 OUTPUTS: 
  g                [double]  n*1 vector, the gap between simulated and sample moments.
  flag             [intger]  empty matrix.

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

0001 function [g,grad,hess,flag] = moment_function(xparams,sample_moments,dataset,options,parallel)
0002 % Evaluates the moment function of the Simulated Moments Method (discrepancy between sample and
0003 % ).
0004 %
0005 % INPUTS:
0006 %  xparams          [double]  p*1 vector of estimated parameters.
0007 %  sample_moments   [double]  n*1 vector of sample moments (n>=p).
0008 %  options          [      ]  Structure defining options for SMM.
0009 %  parallel         [      ]  Structure defining the parallel mode settings (optional).
0010 %
0011 % OUTPUTS:
0012 %  g                [double]  n*1 vector, the gap between simulated and sample moments.
0013 %  flag             [intger]  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 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_
0036 persistent mainStream mainState
0037 persistent priorObjectiveValue
0038 
0039 flag = 1;
0040 grad=[];
0041 hess=[];
0042 
0043 if nargin<5
0044     if isempty(mainStream)
0045         mainStream = RandStream.getDefaultStream;
0046         mainState  = mainStream.State;
0047     else
0048         mainStream.State = mainState;
0049     end
0050 end
0051 
0052 penalty = 0;
0053 for i=1:options.estimated_parameters.nb
0054     if ~isnan(options.estimated_parameters.upper_bound(i)) && xparams(i)>options.estimated_parameters.upper_bound(i)
0055         penalty = penalty + (xparams(i)-options.estimated_parameters.upper_bound(i))^2;
0056     end
0057     if ~isnan(options.estimated_parameters.lower_bound(i)) && xparams(i)<options.estimated_parameters.lower_bound(i)
0058         penalty = penalty + (xparams(i)-options.estimated_parameters.lower_bound(i))^2;
0059     end
0060 end
0061 
0062 if penalty>0
0063     flag = 0;
0064     return;
0065 end
0066 
0067 save('estimated_parameters.mat','xparams');
0068 
0069 % Check for local determinacy of the deterministic steady state.
0070 noprint = options_.noprint; options_.noprint = 1;
0071 [local_determinacy_and_stability,info] = check; options_.noprint = noprint;
0072 if ~local_determinacy_and_stability
0073     flag = 0;
0074     return
0075 end
0076 
0077 simulated_moments = zeros(size(sample_moments));
0078 
0079 % Just to be sure that things don't mess up with persistent variables...
0080 clear perfect_foresight_simulation;
0081 
0082 if nargin<5
0083     for s = 1:options.number_of_simulated_sample
0084         time_series = extended_path([],options.simulated_sample_size,1);
0085         data = time_series(dataset.observed_variables_idx,options.burn_in_periods+1:options.simulated_sample_size);
0086         eval(['tmp = ' options.moments_file_name '(data);'])
0087         simulated_moments = simulated_moments + tmp;
0088         simulated_moments = simulated_moments / options.number_of_simulated_sample;
0089     end
0090 else% parallel mode.
0091     if ~isunix
0092         error('The parallel version of SMM estimation is not implemented for non unix platforms!')
0093     end
0094     job_number = 1;% Remark. First job is for the master.
0095     [Junk,hostname] = unix('hostname --fqdn');
0096     hostname = deblank(hostname);
0097     for i=1:length(parallel)
0098         machine = deblank(parallel(i).machine);
0099         if ~strcmpi(hostname,machine)
0100             % For the slaves on a remote computer.
0101             unix(['scp estimated_parameters.mat ' , parallel(i).login , '@' , machine , ':' parallel(i).folder ' > /dev/null']);
0102         else
0103             if ~strcmpi(pwd,parallel(i).folder)
0104                 % For the slaves on this computer but not in the same directory as the master.
0105                 unix(['cp estimated_parameters.mat ' , parallel(i).folder]);
0106             end
0107         end
0108         for j=1:parallel(i).number_of_jobs
0109             if (strcmpi(hostname,machine) && j>1) || ~strcmpi(hostname,machine)  
0110                 job_number = job_number + 1;
0111                 unix(['ssh -A ' parallel(i).login '@' machine ' ./call_matlab_session.sh job' int2str(job_number) '.m &']);
0112             end
0113         end
0114     end
0115     % Finally the Master do its job
0116     tStartMasterJob = clock;
0117     eval('job1;')
0118     tElapsedMasterJob = etime(clock, tStartMasterJob);
0119     TimeLimit = tElapsedMasterJob*1.2;
0120     % Master waits for the  slaves' output...
0121     tStart = clock;
0122     tElapsed = 0;
0123     while tElapsed<TimeLimit
0124         if ( length(dir('./intermediary_results_from_master_and_slaves/simulated_moments_slave_*.dat'))==job_number )
0125             break
0126         end
0127         tElapsed = etime(clock, tStart);
0128     end
0129     try
0130         tmp = zeros(length(sample_moments),1);
0131         for i=1:job_number
0132             simulated_moments = load(['./intermediary_results_from_master_and_slaves/simulated_moments_slave_' int2str(i) '.dat'],'-ascii');
0133             tmp = tmp + simulated_moments;
0134         end
0135         simulated_moments = tmp / job_number;        
0136     catch
0137         flag = 0;
0138         return
0139     end
0140 end
0141 
0142 g = simulated_moments-sample_moments;