simulated_moments_estimation

 Performs estimation by Simulated Moments Method.

 INPUTS:
  xparam          [double]  p*1 vector of initial values for the estimated parameters. 
  dataset         [      ]  Structure describing the data set.
  options         [      ]  Structure defining options for SMM.
  parallel        [      ]  Structure defining the parallel mode settings (optional). 

 OUTPUTS: 
  param           [double]  p*1 vector of point estimates for the parameters.
  sigma           [double]  p*p covariance matrix of the SMM estimates.

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

0001 function [param,sigma] = simulated_moments_estimation(dataset,options,parallel)
0002 % Performs estimation by Simulated Moments Method.
0003 %
0004 % INPUTS:
0005 %  xparam          [double]  p*1 vector of initial values for the estimated parameters.
0006 %  dataset         [      ]  Structure describing the data set.
0007 %  options         [      ]  Structure defining options for SMM.
0008 %  parallel        [      ]  Structure defining the parallel mode settings (optional).
0009 %
0010 % OUTPUTS:
0011 %  param           [double]  p*1 vector of point estimates for the parameters.
0012 %  sigma           [double]  p*p covariance matrix of the SMM estimates.
0013 %
0014 % SPECIAL REQUIREMENTS
0015 %  The user has to provide a file where the moment conditions are defined.
0016 
0017 % Copyright (C) 2010-2012 Dynare Team
0018 %
0019 % This file is part of Dynare.
0020 %
0021 % Dynare is free software: you can redistribute it and/or modify
0022 % it under the terms of the GNU General Public License as published by
0023 % the Free Software Foundation, either version 3 of the License, or
0024 % (at your option) any later version.
0025 %
0026 % Dynare is distributed in the hope that it will be useful,
0027 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0028 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0029 % GNU General Public License for more details.
0030 %
0031 % You should have received a copy of the GNU General Public License
0032 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0033 
0034 global M_ options_ oo_ estim_params_
0035 
0036 % Load the dataset.
0037 eval(dataset.name);
0038 dataset.data = [];
0039 for v = 1:dataset.number_of_observed_variables
0040     eval(['dataset.data = [ dataset.data , ' dataset.variables(v,:) ' ];'])
0041 end
0042 data = dataset.data(dataset.first_observation:dataset.first_observation+dataset.number_of_observations,:);
0043 
0044 % Compute sample moments and the weighting matrix.
0045 eval(['[sample_moments,long_run_covariance] = ' M_.fname '_moments;'])
0046 weighting_matrix = inv(long_run_covariance);
0047 
0048 % Initialize output.
0049 sigma = [];
0050 param = [];
0051 
0052 % Set options and initial condition.
0053 options.estimated_parameters.list = [];
0054 xparam = [];
0055 if ~isempty(estim_params_.var_exo)
0056     options.estimated_variances.idx = estim_params_.var_exo(:,1);
0057     options.estimated_parameters.list = char(M_.exo_names(options.estimated_variances.idx,:));
0058     options.estimated_parameters.nv = rows(estim_params_.var_exo);
0059     xparam = [xparam; estim_params_.var_exo(:,2)];
0060 end
0061 if ~isempty(estim_params_.param_vals)
0062     options.estimated_parameters.idx = estim_params_.param_vals(:,1);
0063     if isempty(options.estimated_parameters.list)
0064         options.estimated_parameters.list = char(M_.param_names(options.estimated_parameters.idx,:));
0065     else
0066         options.estimated_parameters.list = char(options.estimated_parameters.list,...
0067                                                  M_.param_names(options.estimated_parameters.idx,:));
0068     end
0069     options.estimated_parameters.np = rows(estim_params_.param_vals);
0070     xparam = [xparam; estim_params_.param_vals(:,2)];
0071 end
0072 
0073 options.estimated_parameters.nb = rows(options.estimated_parameters.list);
0074 
0075 options.estimated_parameters.lower_bound = NaN(options.estimated_parameters.nb,1);
0076 options.estimated_parameters.upper_bound = NaN(options.estimated_parameters.nb,1);
0077 
0078 
0079 options.estimated_parameters.lower_bound = [];
0080 options.estimated_parameters.lower_bound = [options.estimated_parameters.lower_bound; ...
0081                     estim_params_.var_exo(:,3); ...
0082                     estim_params_.param_vals(:,3) ];
0083 options.estimated_parameters.upper_bound = [];
0084 options.estimated_parameters.upper_bound = [options.estimated_parameters.upper_bound; ...
0085                     estim_params_.var_exo(:,4); ...
0086                     estim_params_.param_vals(:,4) ];
0087 
0088 options.number_of_simulated_sample = 0;
0089 for i=1:length(parallel)
0090     options.number_of_simulated_sample = options.number_of_simulated_sample + parallel(i).number_of_jobs*parallel(i).number_of_simulations; 
0091 end
0092 
0093 options.observed_variables_idx = dataset.observed_variables_idx;
0094 
0095 % Set up parallel mode if needed.
0096 if nargin>2
0097     if ~isunix
0098         error('The parallel version of SMM estimation is not implemented for non unix platforms!')
0099     end
0100     [junk,hostname] = unix('hostname --fqdn');    
0101     hostname = deblank(hostname);
0102     master_is_running_a_job = 0;
0103     for i=1:length(parallel)
0104         if strcmpi(hostname,parallel(i).machine)
0105             master_is_running_a_job = 1;
0106             break
0107         end
0108     end
0109     if ~master_is_running_a_job
0110         error('Master has to run one job!');
0111     end
0112     if options.optimization_routine>0
0113         estimated_parameters_optimization_path = [NaN;xparam];
0114         save('optimization_path.mat','estimated_parameters_optimization_path');
0115     end
0116     disp(' ')
0117     disp('Master talks to its slaves...')
0118     disp(' ')
0119     % Save the workspace.
0120     save('master_variables.mat','options_','M_','oo_');
0121     % Send the workspace to each remote computer.
0122     disp('')
0123     for i = 1:length(parallel)
0124         if ~strcmpi(hostname,parallel(i).machine)
0125             unix(['scp master_variables.mat ' , parallel(i).login , '@' , parallel(i).machine , ':' parallel(i).folder]);
0126         end
0127     end
0128     for i=1:length(parallel)
0129         % Write a bash script file to execute matlab scripts in the background
0130         fid = fopen('call_matlab_session.sh','w');
0131         fprintf(fid,'#!/bin/sh\n');
0132         fprintf(fid,'unset DISPLAY\n');
0133         fprintf(fid,['cd ' parallel(i).folder '\n']);
0134         fprintf(fid,['nohup ' parallel(i).matlab '/matlab -nodesktop -nodisplay -nojvm < $1 > /dev/null 2>&1\n']);
0135         fprintf(fid,'exit');
0136         fclose(fid);
0137         % Set the permission for this file (has to be executable)
0138         %fileattrib('call_matlab_session.sh','+x','u');
0139         unix(['chmod u+x call_matlab_session.sh']);
0140         % Send the script file on each remote computer
0141         if ~strcmpi(hostname,parallel(i).machine)
0142             unix(['scp call_matlab_session.sh ' , parallel(i).login , '@' , parallel(i).machine , ':~/' ]);
0143         else
0144             unix(['cp call_matlab_session.sh ~/call_matlab_session.sh']);
0145         end
0146     end
0147     % Send the files to each remote computer.
0148     for i = 1:length(parallel)
0149         if ~strcmpi(hostname,parallel(i).machine)
0150             unix(['scp ' M_.fname '_steadystate.m ' , parallel(i).login , '@' , parallel(i).machine , ':' parallel(i).folder]);
0151             unix(['scp ' M_.fname '_moments.m ' , parallel(i).login , '@' , parallel(i).machine , ':' parallel(i).folder]);
0152             unix(['scp ' M_.fname '_static.m ' , parallel(i).login , '@' , parallel(i).machine , ':' parallel(i).folder]);
0153             if exist([M_.fname '_dynamic.c'])
0154                 use_dll_flag = 1;
0155                 unix(['scp ' M_.fname '_dynamic.c ' , parallel(i).login , '@' , parallel(i).machine , ':' parallel(i).folder]);
0156             else
0157                 use_dll_flag = 0;
0158                 unix(['scp ' M_.fname '_dynamic.m ' , parallel(i).login , '@' , parallel(i).machine , ':' parallel(i).folder]);
0159             end
0160         end
0161     end
0162     % If needed, compile dynamic model mex file on each remote computer
0163     if ~strcmpi(hostname,parallel(i).machine) && use_dll_flag
0164         % Write a matlab script that will trigger the compilation of the mex file.
0165         fid = fopen('compile_model.m', 'w');
0166         fprintf(fid,[' eval(''mex -O LDFLAGS=''''-pthread -shared -Wl,--no-undefined'''' ' M_.fname '_dynamic.c'')  ']);
0167         fprintf(fid, '\n exit');
0168         fclose(fid);
0169         for i = 1:length(parallel)
0170             if ~strcmpi(hostname,parallel(i).machine)
0171                 % Send the generated matlab script to the remote computer.
0172                 unix(['scp compile_model.m ' , parallel(i).login , '@' , parallel(i).machine , ':' parallel(i).folder]);
0173                 % Compile the mex file on the remote computer.
0174                 unix(['ssh ' parallel(i).login , '@' , parallel(i).machine ' ./call_matlab_session.sh  compile_model.m']);
0175             end
0176         end
0177     end
0178     % Write the matlab script files for the evaluation of the simulated moment conditions
0179     job = 0;
0180     for i=1:length(parallel)
0181         for j=1:parallel(i).number_of_jobs
0182             job = job+1;
0183             % Create random number streams
0184             write_job(hostname, parallel(i).machine, parallel(i).dynare, ...
0185                       options.simulated_sample_size, length(sample_moments), ...
0186                       dataset.observed_variables_idx, options.estimated_variances.idx', options.estimated_parameters.idx', options.burn_in_periods, [M_.fname '_moments'], parallel(i).number_of_simulations, ...
0187                       parallel(i).number_of_threads_per_job, job, j, options.estimated_parameters.nb, options.estimated_parameters.nv, ...
0188                       options.estimated_parameters.np);
0189             if ~strcmpi(hostname,parallel(i).machine)
0190                 unix(['scp ' , 'job' , int2str(job) , '.m ' , parallel(i).login , '@' , parallel(i).machine , ':' parallel(i).folder ]);
0191             end
0192         end
0193     end
0194     disp(' ')
0195     disp('... And slaves do as ordered.')
0196     disp(' ')
0197     if exist('intermediary_results_from_master_and_slaves','dir')
0198         unix('rm -rf intermediary_results_from_master_and_slaves');
0199     end
0200     unix('mkdir intermediary_results_from_master_and_slaves');
0201     unix('chmod -R u+x intermediary_results_from_master_and_slaves');
0202 end
0203 
0204 disp('');
0205 
0206 if options.optimization_routine==1
0207     % Set options for csminwel.
0208     H0 = 1e-4*eye(options.estimated_parameters.nb);
0209     ct = 1e-4;
0210     it = 1000;
0211     vb = 2;
0212     % Minimization of the objective function.
0213     if nargin==2
0214         [fval,param,grad,hessian_csminwel,itct,fcount,retcodehat] = ...
0215             csminwel1('smm_objective',xparam,H0,[],ct,it,2,options_.gradient_epsilon,sample_moments,weighting_matrix,options);    
0216     elseif nargin>2
0217         [fval,param,grad,hessian_csminwel,itct,fcount,retcodehat] = ...
0218             csminwel1('smm_objective',xparam,H0,[],ct,it,2,options_.gradient_epsilon,sample_moments,weighting_matrix,options,parallel);
0219     end
0220 elseif options.optimization_routine==2
0221     optim_options = optimset('display','iter','MaxFunEvals',1000000,'MaxIter',6000,'TolFun',1e-4,'TolX',1e-4);
0222     if isfield(options_,'optim_opt')
0223         eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
0224     end
0225     if nargin==2
0226         [param,fval,exitflag] = fminsearch('smm_objective',xparam,optim_options,sample_moments,weighting_matrix,options);
0227     else
0228         [param,fval,exitflag] = fminsearch('smm_objective',xparam,optim_options,sample_moments,weighting_matrix,options,parallel);
0229     end
0230 elseif options.optimization_routine==0% Compute the variance of the SMM estimator
0231     load('optimization_path.mat');
0232     tmp = sortrows(estimated_parameters_optimization_path',1);
0233     param = tmp(1,2:end)';
0234     % Compute gradient of the moment function (distance between sample and simulated moments).
0235     [F,G] = dynare_gradient('moment_function',param,options_.gradient_epsilon,sample_moments,dataset,options,parallel);
0236     V = (1+1/options.number_of_simulated_sample)*G'*long_run_covariance*G;
0237     [param,diag(V)]
0238 elseif options.optimization_routine<0
0239     T = -options.optimization_routine;% length of the simulated time series.
0240     time_series = extended_path(oo_.steady_state,T,1);
0241     save time_series.mat;
0242 end
0243 
0244 
0245 function write_job(hostname, remotename, dynare_path, sample_size, number_of_moments, observed_variables_idx, variance_idx, parameters_idx, burn_in_periods, moments_file_name, number_of_simulations,threads_per_job, slave_number, job_number,nb,nv,np)
0246 
0247 fid = fopen(['job' int2str(slave_number) '.m'],'w');
0248 
0249 fprintf(fid,['% Generated by ' hostname '.\n\n']);
0250 
0251 if ( strcmpi(hostname,remotename) && (job_number>1) )  || ~strcmpi(hostname,remotename) 
0252     fprintf(fid,'load(''master_variables'');\n');
0253     fprintf(fid,'assignin(''base'',''M_'',M_);\n');
0254     fprintf(fid,'assignin(''base'',''oo_'',oo_);\n');
0255     fprintf(fid,'assignin(''base'',''options_'',options_);\n\n');
0256 end
0257 
0258 if ( strcmpi(hostname,remotename) && (job_number>1) ) || ~strcmpi(hostname,remotename)
0259     fprintf(fid,['addpath ' dynare_path '\n']);
0260     fprintf(fid,['dynare_config;\n\n']);
0261 end
0262 
0263 fprintf(fid,['simulated_moments = zeros(' int2str(number_of_moments) ',1);\n\n']);
0264 
0265 fprintf(fid,'load(''estimated_parameters.mat'');\n');
0266 fprintf(fid,['M_.params([' num2str(parameters_idx)  ']) = xparams(' int2str(nv) '+1:' int2str(nb) ');\n\n']);
0267 fprintf(fid,'tmp = diag(M_.Sigma_e);')
0268 fprintf(fid,['tmp([' num2str(variance_idx)  ']) = xparams(1:' int2str(nv) ').^2;\n\n']);
0269 fprintf(fid,'M_.Sigma_e = diag(tmp);')
0270 
0271 fprintf(fid,['stream=RandStream(''mt19937ar'',''Seed'',' int2str(slave_number) ');\n']);
0272 if matlab_ver_less_than('7.12')
0273       fprintf(fid,['RandStream.setDefaultStream(stream);\n\n']);
0274 else
0275       fprintf(fid,['RandStream.setGlobalStream(stream);\n\n']);
0276 end
0277 
0278 fprintf(fid,['maxNumCompThreads(' int2str(threads_per_job) ');\n\n']);
0279 
0280 fprintf(fid,['for s = 1:' int2str(number_of_simulations) '\n'] );
0281 fprintf(fid,['    time_series = extended_path([],' int2str(sample_size) ',1);\n']);
0282 fprintf(fid,['    data = time_series([' int2str(observed_variables_idx) '],' int2str(burn_in_periods) '+1:' int2str(sample_size) ');\n']);
0283 fprintf(fid,['    eval(''tmp = ' moments_file_name '(data);'');\n']);
0284 fprintf(fid,['    simulated_moments = simulated_moments + tmp;\n']);
0285 fprintf(fid,['end;\n\n']);
0286 
0287 fprintf(fid,['simulated_moments = simulated_moments/' int2str(number_of_simulations) ';\n']);
0288 fprintf(fid,['save(''simulated_moments_slave_' int2str(slave_number) '.dat'',''simulated_moments'',''-ascii'');\n']);
0289 
0290 if ~strcmpi(hostname,remotename)
0291     fprintf(fid,['unix(''scp simulated_moments_slave_' int2str(slave_number) '.dat ' hostname ':' pwd '/intermediary_results_from_master_and_slaves '');\n']);
0292     fprintf(fid,['unix(''rm simulated_moments_slave_' int2str(slave_number) '.dat'');\n']);
0293 else
0294     fprintf(fid,['unix(''cp simulated_moments_slave_' int2str(slave_number) '.dat '  'intermediary_results_from_master_and_slaves '');\n']);
0295     fprintf(fid,['unix(''rm simulated_moments_slave_' int2str(slave_number) '.dat'');\n']);
0296 end
0297 
0298 if ((job_number>1) && strcmpi(hostname,remotename)) || ~strcmpi(hostname,remotename) 
0299     fprintf(fid,'exit');
0300 end
0301 
0302 fclose(fid);