prior_posterior_statistics_core

 PARALLEL CONTEXT
 Core functionality for prior_posterior.m function, which can be parallelized.
 See also the comment in random_walk_metropolis_hastings_core.m funtion.

 INPUTS
   See See the comment in random_walk_metropolis_hastings_core.m funtion.

 OUTPUTS
 o myoutput  [struc]
  Contained OutputFileName_smooth;
                          _update;
                          _inno;
                          _error;
                          _filter_step_ahead;
                          _param;
                          _forc_mean;
                          _forc_point

 ALGORITHM
   Portion of prior_posterior.m function.
 This file is part of Dynare.

 SPECIAL REQUIREMENTS.
   None.

0001 function myoutput=prior_posterior_statistics_core(myinputs,fpar,B,whoiam, ThisMatlab)
0002 % PARALLEL CONTEXT
0003 % Core functionality for prior_posterior.m function, which can be parallelized.
0004 % See also the comment in random_walk_metropolis_hastings_core.m funtion.
0005 %
0006 % INPUTS
0007 %   See See the comment in random_walk_metropolis_hastings_core.m funtion.
0008 %
0009 % OUTPUTS
0010 % o myoutput  [struc]
0011 %  Contained OutputFileName_smooth;
0012 %                          _update;
0013 %                          _inno;
0014 %                          _error;
0015 %                          _filter_step_ahead;
0016 %                          _param;
0017 %                          _forc_mean;
0018 %                          _forc_point
0019 %
0020 % ALGORITHM
0021 %   Portion of prior_posterior.m function.
0022 % This file is part of Dynare.
0023 %
0024 % SPECIAL REQUIREMENTS.
0025 %   None.
0026 
0027 % Copyright (C) 2005-2011 Dynare Team
0028 %
0029 % This file is part of Dynare.
0030 %
0031 % Dynare is free software: you can redistribute it and/or modify
0032 % it under the terms of the GNU General Public License as published by
0033 % the Free Software Foundation, either version 3 of the License, or
0034 % (at your option) any later version.
0035 %
0036 % Dynare is distributed in the hope that it will be useful,
0037 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0038 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0039 % GNU General Public License for more details.
0040 %
0041 % You should have received a copy of the GNU General Public License
0042 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0043 
0044 global options_ oo_ M_ bayestopt_ estim_params_
0045 
0046 if nargin<4,
0047     whoiam=0;
0048 end
0049 
0050 global options_ oo_ M_ bayestopt_ estim_params_
0051 
0052 
0053 % Reshape 'myinputs' for local computation.
0054 % In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
0055 
0056 type=myinputs.type;
0057 run_smoother=myinputs.run_smoother;
0058 gend=myinputs.gend;
0059 Y=myinputs.Y;
0060 data_index=myinputs.data_index;
0061 missing_value=myinputs.missing_value;
0062 varobs=myinputs.varobs;
0063 irun=myinputs.irun;
0064 endo_nbr=myinputs.endo_nbr;
0065 nvn=myinputs.nvn;
0066 naK=myinputs.naK;
0067 horizon=myinputs.horizon;
0068 iendo=myinputs.iendo;
0069 if horizon
0070     i_last_obs=myinputs.i_last_obs;
0071     IdObs=myinputs.IdObs;
0072     MAX_nforc1=myinputs.MAX_nforc1;
0073     MAX_nforc2=myinputs.MAX_nforc2;
0074 end
0075 if naK
0076     MAX_naK=myinputs.MAX_naK;
0077 end
0078 
0079 exo_nbr=myinputs.exo_nbr;
0080 maxlag=myinputs.maxlag;
0081 MAX_nsmoo=myinputs.MAX_nsmoo;
0082 MAX_ninno=myinputs.MAX_ninno;
0083 MAX_nerro = myinputs.MAX_nerro;
0084 MAX_nruns=myinputs.MAX_nruns;
0085 MAX_momentsno = myinputs.MAX_momentsno;
0086 ifil=myinputs.ifil;
0087 
0088 if ~strcmpi(type,'prior'),
0089     x=myinputs.x;
0090     if strcmpi(type,'posterior'),
0091         logpost=myinputs.logpost;
0092     end
0093 end
0094 if whoiam
0095     Parallel=myinputs.Parallel;
0096 end
0097 
0098 % DirectoryName = myinputs.DirectoryName;
0099 if strcmpi(type,'posterior')
0100     DirectoryName = CheckPath('metropolis',M_.dname);
0101 elseif strcmpi(type,'gsa')
0102     if options_.opt_gsa.pprior
0103         DirectoryName = CheckPath(['gsa',filesep,'prior'],M_.dname);
0104     else
0105         DirectoryName = CheckPath(['gsa',filesep,'mc'],M_.dname);
0106     end
0107 elseif strcmpi(type,'prior')
0108     DirectoryName = CheckPath('prior',M_.dname);
0109 end
0110 
0111 RemoteFlag = 0;
0112 if whoiam
0113     if Parallel(ThisMatlab).Local==0,
0114         RemoteFlag =1;
0115     end
0116     ifil=ifil(:,whoiam);
0117     prct0={0,whoiam,Parallel(ThisMatlab)};
0118 else
0119     prct0=0;
0120 end
0121 h = dyn_waitbar(prct0,['Taking ',type,' subdraws...']);
0122 
0123 if RemoteFlag==1,
0124     OutputFileName_smooth = {};
0125     OutputFileName_update = {};
0126     OutputFileName_inno = {};
0127     OutputFileName_error = {};
0128     OutputFileName_filter_step_ahead = {};
0129     OutputFileName_param = {};
0130     OutputFileName_forc_mean = {};
0131     OutputFileName_forc_point = {};
0132     % OutputFileName_moments = {};
0133 end
0134 
0135 for b=fpar:B
0136 
0137     %    [deep, logpo] = GetOneDraw(type);
0138     %    set_all_parameters(deep);
0139     %    dr = resol(oo_.steady_state,0);
0140     if strcmpi(type,'prior')
0141 
0142         [deep, logpo] = GetOneDraw(type);
0143 
0144     else
0145         deep = x(b,:);
0146         if strcmpi(type,'posterior')
0147             logpo = logpost(b);
0148         else
0149             logpo = evaluate_posterior_kernel(deep');
0150         end
0151     end
0152     set_all_parameters(deep);
0153     [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
0154 
0155     if run_smoother
0156         [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK] = ...
0157             DsgeSmoother(deep,gend,Y,data_index,missing_value);
0158 
0159         if options_.loglinear
0160             stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ...
0161                 repmat(log(SteadyState(dr.order_var)),1,gend);
0162             stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ...
0163                 repmat(log(SteadyState(dr.order_var)),1,gend);
0164         else
0165             stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ...
0166                 repmat(SteadyState(dr.order_var),1,gend);
0167             stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ...
0168                 repmat(SteadyState(dr.order_var),1,gend);
0169         end
0170         stock_innov(:,:,irun(2))  = etahat;
0171         if nvn
0172             stock_error(:,:,irun(3))  = epsilonhat;
0173         end
0174         if naK
0175             stock_filter_step_ahead(:,dr.order_var,:,irun(4)) = aK(options_.filter_step_ahead,1:endo_nbr,:);
0176         end
0177 
0178         if horizon
0179             yyyy = alphahat(iendo,i_last_obs);
0180             yf = forcst2a(yyyy,dr,zeros(horizon,exo_nbr));
0181             if options_.prefilter == 1
0182                 yf(:,IdObs) = yf(:,IdObs)+repmat(bayestopt_.mean_varobs', ...
0183                                                  horizon+maxlag,1);
0184             end
0185             yf(:,IdObs) = yf(:,IdObs)+(gend+[1-maxlag:horizon]')*trend_coeff';
0186             if options_.loglinear == 1
0187                 yf = yf+repmat(log(SteadyState'),horizon+maxlag,1);
0188             else
0189                 yf = yf+repmat(SteadyState',horizon+maxlag,1);
0190             end
0191             yf1 = forcst2(yyyy,horizon,dr,1);
0192             if options_.prefilter == 1
0193                 yf1(:,IdObs,:) = yf1(:,IdObs,:)+ ...
0194                     repmat(bayestopt_.mean_varobs',[horizon+maxlag,1,1]);
0195             end
0196             yf1(:,IdObs,:) = yf1(:,IdObs,:)+repmat((gend+[1-maxlag:horizon]')* ...
0197                                                    trend_coeff',[1,1,1]);
0198             if options_.loglinear == 1
0199                 yf1 = yf1 + repmat(log(SteadyState'),[horizon+maxlag,1,1]);
0200             else
0201                 yf1 = yf1 + repmat(SteadyState',[horizon+maxlag,1,1]);
0202             end
0203 
0204             stock_forcst_mean(:,:,irun(6)) = yf';
0205             stock_forcst_point(:,:,irun(7)) = yf1';
0206         end
0207 
0208     end
0209     stock_param(irun(5),:) = deep;
0210     stock_logpo(irun(5),1) = logpo;
0211     stock_ys(irun(5),:) = SteadyState';
0212 
0213     irun = irun +  ones(7,1);
0214 
0215 
0216     if irun(1) > MAX_nsmoo || b == B
0217         stock = stock_smooth(:,:,1:irun(1)-1);
0218         ifil(1) = ifil(1) + 1;
0219         save([DirectoryName '/' M_.fname '_smooth' int2str(ifil(1)) '.mat'],'stock');
0220 
0221         stock = stock_update(:,:,1:irun(1)-1);
0222         save([DirectoryName '/' M_.fname '_update' int2str(ifil(1)) '.mat'],'stock');
0223         if RemoteFlag==1,
0224             OutputFileName_smooth = [OutputFileName_smooth; {[DirectoryName filesep], [M_.fname '_smooth' int2str(ifil(1)) '.mat']}];
0225             OutputFileName_update = [OutputFileName_update; {[DirectoryName filesep], [M_.fname '_update' int2str(ifil(1)) '.mat']}];
0226         end
0227         irun(1) = 1;
0228     end
0229 
0230     if irun(2) > MAX_ninno || b == B
0231         stock = stock_innov(:,:,1:irun(2)-1);
0232         ifil(2) = ifil(2) + 1;
0233         save([DirectoryName '/' M_.fname '_inno' int2str(ifil(2)) '.mat'],'stock');
0234         if RemoteFlag==1,
0235             OutputFileName_inno = [OutputFileName_inno; {[DirectoryName filesep], [M_.fname '_inno' int2str(ifil(2)) '.mat']}];
0236         end
0237         irun(2) = 1;
0238     end
0239 
0240     if nvn && (irun(3) > MAX_nerro || b == B)
0241         stock = stock_error(:,:,1:irun(3)-1);
0242         ifil(3) = ifil(3) + 1;
0243         save([DirectoryName '/' M_.fname '_error' int2str(ifil(3)) '.mat'],'stock');
0244         if RemoteFlag==1,
0245             OutputFileName_error = [OutputFileName_error; {[DirectoryName filesep], [M_.fname '_error' int2str(ifil(3)) '.mat']}];
0246         end
0247         irun(3) = 1;
0248     end
0249 
0250     if naK && (irun(4) > MAX_naK || b == B)
0251         stock = stock_filter_step_ahead(:,:,:,1:irun(4)-1);
0252         ifil(4) = ifil(4) + 1;
0253         save([DirectoryName '/' M_.fname '_filter_step_ahead' int2str(ifil(4)) '.mat'],'stock');
0254         if RemoteFlag==1,
0255             OutputFileName_filter_step_ahead = [OutputFileName_filter_step_ahead; {[DirectoryName filesep], [M_.fname '_filter_step_ahead' int2str(ifil(4)) '.mat']}];
0256         end
0257         irun(4) = 1;
0258     end
0259 
0260     if irun(5) > MAX_nruns || b == B
0261         stock = stock_param(1:irun(5)-1,:);
0262         ifil(5) = ifil(5) + 1;
0263         save([DirectoryName '/' M_.fname '_param' int2str(ifil(5)) '.mat'],'stock','stock_logpo','stock_ys');
0264         if RemoteFlag==1,
0265             OutputFileName_param = [OutputFileName_param; {[DirectoryName filesep], [M_.fname '_param' int2str(ifil(5)) '.mat']}];
0266         end
0267         irun(5) = 1;
0268     end
0269 
0270     if horizon && (irun(6) > MAX_nforc1 || b == B)
0271         stock = stock_forcst_mean(:,:,1:irun(6)-1);
0272         ifil(6) = ifil(6) + 1;
0273         save([DirectoryName '/' M_.fname '_forc_mean' int2str(ifil(6)) '.mat'],'stock');
0274         if RemoteFlag==1,
0275             OutputFileName_forc_mean = [OutputFileName_forc_mean; {[DirectoryName filesep], [M_.fname '_forc_mean' int2str(ifil(6)) '.mat']}];
0276         end
0277         irun(6) = 1;
0278     end
0279 
0280     if horizon && (irun(7) > MAX_nforc2 ||  b == B)
0281         stock = stock_forcst_point(:,:,1:irun(7)-1);
0282         ifil(7) = ifil(7) + 1;
0283         save([DirectoryName '/' M_.fname '_forc_point' int2str(ifil(7)) '.mat'],'stock');
0284         if RemoteFlag==1,
0285             OutputFileName_forc_point = [OutputFileName_forc_point; {[DirectoryName filesep], [M_.fname '_forc_point' int2str(ifil(7)) '.mat']}];
0286         end
0287         irun(7) = 1;
0288     end
0289 
0290     % if moments_varendo && (irun(8) > MAX_momentsno || b == B)
0291     %    stock = stock_moments(1:irun(8)-1);
0292     %    ifil(8) = ifil(8) + 1;
0293     %    save([DirectoryName '/' M_.fname '_moments' int2str(ifil(8)) '.mat'],'stock');
0294     %    if RemoteFlag==1,
0295     %    OutputFileName_moments = [OutputFileName_moments; {[DirectoryName filesep], [M_.fname '_moments' int2str(ifil(8)) '.mat']}];
0296     %    end
0297     %    irun(8) = 1;
0298     % end
0299 
0300     %   DirectoryName=TempPath;
0301 
0302 
0303     dyn_waitbar((b-fpar+1)/(B-fpar+1),h);
0304 
0305 end
0306 
0307 myoutput.ifil=ifil;
0308 if RemoteFlag==1,
0309     myoutput.OutputFileName = [OutputFileName_smooth;
0310                         OutputFileName_update;
0311                         OutputFileName_inno;
0312                         OutputFileName_error;
0313                         OutputFileName_filter_step_ahead;
0314                         OutputFileName_param;
0315                         OutputFileName_forc_mean;
0316                         OutputFileName_forc_point];
0317     % OutputFileName_moments];
0318 end
0319 
0320 dyn_waitbar_close(h);
0321