From 2f717b5adc5a87f663c5c080f2963d1f65d1933e Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Thu, 15 Dec 2016 10:36:16 +0100
Subject: [PATCH] Eliminate global variables from shock_decomposition.m

---
 matlab/DsgeSmoother.m                    | 13 ++++++++++---
 matlab/dynare_estimation_1.m             |  4 ++--
 matlab/evaluate_smoother.m               | 14 ++++++++------
 matlab/gsa/filt_mc_.m                    |  4 ++--
 matlab/imcforecast.m                     |  4 ++--
 matlab/prior_posterior_statistics_core.m |  4 ++--
 matlab/shock_decomposition.m             |  6 ++++--
 preprocessor/ComputingTasks.cc           |  4 ++--
 8 files changed, 32 insertions(+), 21 deletions(-)

diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 6efd3a455..4dd5711af 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -1,4 +1,4 @@
-function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition,state_uncertainty] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value)
+function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition,state_uncertainty,M_,oo_,options_,bayestopt_] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_)
 % Estimation of the smoothed variables and innovations. 
 % 
 % INPUTS 
@@ -7,6 +7,11 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,de
 %   o data          [double]   (n*T) matrix of data.
 %   o data_index    [cell]      1*smpl cell of column vectors of indices.
 %   o missing_value 1 if missing values, 0 otherwise
+%   o M_            [structure] decribing the model
+%   o oo_           [structure] storing the results
+%   o options_      [structure] describing the options
+%   o bayestopt_    [structure] describing the priors
+%   o estim_params_ [structure] characterizing parameters to be estimated
 %  
 % OUTPUTS
 %   o alphahat      [double]  (m*T) matrix, smoothed endogenous variables (a_{t|T})  (decision-rule order)
@@ -27,6 +32,10 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,de
 %   o trend_addition [double] (n*T) pure trend component; stored in options_.varobs order         
 %   o state_uncertainty [double] (K,K,T) array, storing the uncertainty
 %                                   about the smoothed state (decision-rule order)
+%   o M_            [structure] decribing the model
+%   o oo_           [structure] storing the results
+%   o options_      [structure] describing the options
+%   o bayestopt_    [structure] describing the priors
 %  
 % Notes:
 %   m:  number of endogenous variables (M_.endo_nbr)
@@ -66,8 +75,6 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,de
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global bayestopt_ M_ oo_ estim_params_ options_
-
 alphahat        = [];
 etahat  = [];
 epsilonhat      = [];
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index ff9e86836..e3704bea8 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -163,7 +163,7 @@ end
 
 if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.mh_posterior_mode_estimation==0
     if options_.smoother == 1
-        [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value);
+        [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,options_,bayestopt_] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
         [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
     end
     %reset qz_criterium
@@ -483,7 +483,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
                                                       > 0) && options_.load_mh_file)) ...
     || ~options_.smoother ) && options_.partial_information == 0  % to be fixed
     %% ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothes variables
-    [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty] = DsgeSmoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state);
+    [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,options_,bayestopt_] = DsgeSmoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_);
     [oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
 
     if ~options_.nograph,
diff --git a/matlab/evaluate_smoother.m b/matlab/evaluate_smoother.m
index 6e420d95a..72cf57931 100644
--- a/matlab/evaluate_smoother.m
+++ b/matlab/evaluate_smoother.m
@@ -1,11 +1,15 @@
-function [oo_, Smoothed_variables_declaration_order_deviation_form]=evaluate_smoother(parameters,var_list)
+function [oo_, Smoothed_variables_declaration_order_deviation_form]=evaluate_smoother(parameters,var_list,M_,oo_,options_,bayestopt_,estim_params_)
 % Evaluate the smoother at parameters.
 %
 % INPUTS
 %    o parameters  a string ('posterior mode','posterior mean','posterior median','prior mode','prior mean','mle_mode') or a vector of values for
 %                  the (estimated) parameters of the model.
 %    o var_list    subset of endogenous variables
-%
+%    o M_          [structure]  Definition of the model
+%    o oo_         [structure]  Storage of results
+%    o options_    [structure]  Options
+%    o bayestopt_  [structure]  describing the priors
+%    o estim_params_ [structure] characterizing parameters to be estimated
 %
 % OUTPUTS
 %    o oo       [structure]  results:
@@ -47,8 +51,6 @@ function [oo_, Smoothed_variables_declaration_order_deviation_form]=evaluate_smo
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global options_ M_ bayestopt_ oo_ estim_params_   % estim_params_ may be emty
-
 persistent dataset_ dataset_info
 
 %store qz_criterium
@@ -97,8 +99,8 @@ if ischar(parameters)
     end
 end
 
-[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty] = ...
-    DsgeSmoother(parameters,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state);
+[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,options_,bayestopt_] = ...
+    DsgeSmoother(parameters,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_);
 [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
 
 if nargout==2
diff --git a/matlab/gsa/filt_mc_.m b/matlab/gsa/filt_mc_.m
index 2d3050926..3c79676d8 100644
--- a/matlab/gsa/filt_mc_.m
+++ b/matlab/gsa/filt_mc_.m
@@ -163,7 +163,7 @@ if ~loadSA,
             yss(i,:,:)=repmat(sto_ys(:,js(i))',[gend,1]);
         end
         if exist('xparam1','var')
-            [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value);
+            [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
             y0 = transpose( squeeze(aK(1,jxj,1:gend)));% + kron(ys_mode(js),ones(1,gend)));
             yobs = transpose( ahat(jxj,:));% + kron(ys_mode(js),ones(1,gend)));
             rmse_mode = sqrt(mean((yobs(istart:end,:)-y0(istart:end,:)).^2));
@@ -196,7 +196,7 @@ if ~loadSA,
         end
         if exist('xparam1_mean','var')
             %eval(['rmse_pmean(i) = sqrt(mean((',vj,'(fobs-1+istart:fobs-1+nobs)-y0M(istart:end-1)).^2));'])
-            [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1_mean,gend,Y,data_index,missing_value);
+            [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1_mean,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
             y0 = transpose( squeeze(aK(1,jxj,1:gend)));% + kron(ys_mean(js),ones(1,gend)));
             yobs = transpose( ahat(jxj,:));% + kron(ys_mean(js),ones(1,gend)));
             rmse_pmean = sqrt(mean((yobs(istart:end,:)-y0(istart:end,:)).^2));
diff --git a/matlab/imcforecast.m b/matlab/imcforecast.m
index 5a9dbb3e6..f40475e48 100644
--- a/matlab/imcforecast.m
+++ b/matlab/imcforecast.m
@@ -43,7 +43,7 @@ function imcforecast(constrained_paths, constrained_vars, options_cond_fcst)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global options_ oo_ M_ bayestopt_
+global options_ oo_ M_ bayestopt_ estim_params_
 
 if ~isfield(options_cond_fcst,'parameter_set') || isempty(options_cond_fcst.parameter_set)
     if isfield(oo_,'posterior_mode')
@@ -124,7 +124,7 @@ if estimated_model
     qz_criterium_old=options_.qz_criterium;
     options_=select_qz_criterium_value(options_);
     options_smoothed_state_uncertainty_old = options_.smoothed_state_uncertainty;
-    [atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff,aK,T,R,P,PK,decomp,trend_addition] = DsgeSmoother(xparam,gend,data,data_index,missing_value);
+    [atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff,aK,T,R,P,PK,decomp,trend_addition,state_uncertainty,M_,oo_,options_,bayestopt_] = DsgeSmoother(xparam,gend,data,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
     options_.smoothed_state_uncertainty = options_smoothed_state_uncertainty_old;
     %get constant part
     if options_.noconstant
diff --git a/matlab/prior_posterior_statistics_core.m b/matlab/prior_posterior_statistics_core.m
index df250d79a..3e8e31436 100644
--- a/matlab/prior_posterior_statistics_core.m
+++ b/matlab/prior_posterior_statistics_core.m
@@ -204,8 +204,8 @@ for b=fpar:B
 
     if run_smoother
         [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
-        [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,junk1,junk2,P,junk4,junk5,trend_addition,state_uncertainty] = ...
-            DsgeSmoother(deep,gend,Y,data_index,missing_value);
+        [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,junk1,junk2,P,junk4,junk5,trend_addition,state_uncertainty,M_,oo_,options_,bayestopt_] = ...
+            DsgeSmoother(deep,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
         
         stock_trend_coeff(options_.varobs_id,irun(9))=trend_coeff;
         stock_smoothed_trend(IdObs,:,irun(11))=trend_addition;
diff --git a/matlab/shock_decomposition.m b/matlab/shock_decomposition.m
index b56a591a1..91205a037 100644
--- a/matlab/shock_decomposition.m
+++ b/matlab/shock_decomposition.m
@@ -1,4 +1,4 @@
-function oo_ = shock_decomposition(M_,oo_,options_,varlist)
+function oo_ = shock_decomposition(M_,oo_,options_,varlist,bayestopt_,estim_params_)
 % function z = shock_decomposition(M_,oo_,options_,varlist)
 % Computes shocks contribution to a simulated trajectory. The field set is
 % oo_.shock_decomposition. It is a n_var by nshock+2 by nperiods array. The
@@ -12,6 +12,8 @@ function oo_ = shock_decomposition(M_,oo_,options_,varlist)
 %    oo_:         [structure]  Storage of results
 %    options_:    [structure]  Options
 %    varlist:     [char]       List of variables
+%    bayestopt_:  [structure]  describing the priors
+%    estim_params_: [structure] characterizing parameters to be estimated
 %
 % OUTPUTS
 %    oo_:         [structure]  Storage of results
@@ -64,7 +66,7 @@ if isempty(parameter_set)
     end
 end
 
-[oo,Smoothed_Variables_deviation_from_mean] = evaluate_smoother(parameter_set,varlist);
+[oo,Smoothed_Variables_deviation_from_mean] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_);
 
 % reduced form
 dr = oo.dr;
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index 17414e763..1a4426208 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -1642,7 +1642,7 @@ ShockDecompositionStatement::writeOutput(ostream &output, const string &basename
 {
   options_list.writeOutput(output);
   symbol_list.writeOutput("var_list_", output);
-  output << "oo_ = shock_decomposition(M_,oo_,options_,var_list_);" << endl;
+  output << "oo_ = shock_decomposition(M_,oo_,options_,var_list_,bayestopt_,estim_params_);" << endl;
 }
 
 ConditionalForecastStatement::ConditionalForecastStatement(const OptionsList &options_list_arg) :
@@ -3171,7 +3171,7 @@ CalibSmootherStatement::writeOutput(ostream &output, const string &basename, boo
   symbol_list.writeOutput("var_list_", output);
   output << "options_.smoother = 1;" << endl;
   output << "options_.order = 1;" << endl;
-  output << "evaluate_smoother('calibration',var_list_);" << endl;
+  output << "evaluate_smoother('calibration',var_list_,M_,oo_,options_,bayestopt_,estim_params_);" << endl;
 }
 
 ExtendedPathStatement::ExtendedPathStatement(const OptionsList &options_list_arg)
-- 
GitLab