prior_posterior_statistics_core.m 10.9 KB
Newer Older
Marco Ratto's avatar
Marco Ratto committed
1
function myoutput=prior_posterior_statistics_core(myinputs,fpar,B,whoiam, ThisMatlab)
2 3 4 5
% 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.
%
6
% INPUTS
7
%   See See the comment in random_walk_metropolis_hastings_core.m funtion.
8
%
9 10
% OUTPUTS
% o myoutput  [struc]
11 12 13 14 15
%  Contained OutputFileName_smooth;
%                          _update;
%                          _inno;
%                          _error;
%                          _filter_step_ahead;
16 17 18
%                          _param;
%                          _forc_mean;
%                          _forc_point
19
%
20 21
% ALGORITHM
%   Portion of prior_posterior.m function.
22 23
% This file is part of Dynare.
%
24 25
% SPECIAL REQUIREMENTS.
%   None.
26

27
% Copyright (C) 2005-2013 Dynare Team
28 29 30
%
% This file is part of Dynare.
%
31 32 33 34 35 36 37 38 39 40 41 42
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
Marco Ratto's avatar
Marco Ratto committed
43

44 45
global options_ oo_ M_ bayestopt_ estim_params_

Marco Ratto's avatar
Marco Ratto committed
46 47 48 49
if nargin<4,
    whoiam=0;
end

50 51 52 53
% Reshape 'myinputs' for local computation.
% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:

type=myinputs.type;
54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
run_smoother=myinputs.run_smoother;
gend=myinputs.gend;
Y=myinputs.Y;
data_index=myinputs.data_index;
missing_value=myinputs.missing_value;
varobs=myinputs.varobs;
irun=myinputs.irun;
endo_nbr=myinputs.endo_nbr;
nvn=myinputs.nvn;
naK=myinputs.naK;
horizon=myinputs.horizon;
iendo=myinputs.iendo;
if horizon
    i_last_obs=myinputs.i_last_obs;
    IdObs=myinputs.IdObs;
69 70
    MAX_nforc1=myinputs.MAX_nforc1;
    MAX_nforc2=myinputs.MAX_nforc2;
71
end
72 73 74 75
if naK
    MAX_naK=myinputs.MAX_naK;
end

76 77 78 79
exo_nbr=myinputs.exo_nbr;
maxlag=myinputs.maxlag;
MAX_nsmoo=myinputs.MAX_nsmoo;
MAX_ninno=myinputs.MAX_ninno;
80
MAX_nerro = myinputs.MAX_nerro;
81
MAX_nruns=myinputs.MAX_nruns;
82
MAX_momentsno = myinputs.MAX_momentsno;
83
ifil=myinputs.ifil;
Marco Ratto's avatar
Marco Ratto committed
84

85
if ~strcmpi(type,'prior'),
86
    x=myinputs.x;
87 88 89
    if strcmpi(type,'posterior'),
        logpost=myinputs.logpost;
    end
90
end
91 92 93
if whoiam
    Parallel=myinputs.Parallel;
end
Marco Ratto's avatar
Marco Ratto committed
94

95 96
% DirectoryName = myinputs.DirectoryName;
if strcmpi(type,'posterior')
97
    DirectoryName = CheckPath('metropolis',M_.dname);
98 99
elseif strcmpi(type,'gsa')
    if options_.opt_gsa.pprior
100
        DirectoryName = CheckPath(['gsa',filesep,'prior'],M_.dname);
101
    else
102
        DirectoryName = CheckPath(['gsa',filesep,'mc'],M_.dname);
103 104
    end
elseif strcmpi(type,'prior')
105
    DirectoryName = CheckPath('prior',M_.dname);
106
end
Marco Ratto's avatar
Marco Ratto committed
107 108 109

RemoteFlag = 0;
if whoiam
110 111
    if Parallel(ThisMatlab).Local==0,
        RemoteFlag =1;
112
    end
Marco Ratto's avatar
Marco Ratto committed
113
    ifil=ifil(:,whoiam);
114
    prct0={0,whoiam,Parallel(ThisMatlab)};
115
else
116
    prct0=0;
Marco Ratto's avatar
Marco Ratto committed
117
end
118
h = dyn_waitbar(prct0,['Taking ',type,' subdraws...']);
Marco Ratto's avatar
Marco Ratto committed
119 120

if RemoteFlag==1,
121 122 123 124 125 126 127 128 129
    OutputFileName_smooth = {};
    OutputFileName_update = {};
    OutputFileName_inno = {};
    OutputFileName_error = {};
    OutputFileName_filter_step_ahead = {};
    OutputFileName_param = {};
    OutputFileName_forc_mean = {};
    OutputFileName_forc_point = {};
    % OutputFileName_moments = {};
Marco Ratto's avatar
Marco Ratto committed
130 131
end

132 133 134 135
%initialize arrays
if run_smoother
  stock_smooth=NaN(endo_nbr,gend,MAX_nsmoo);
  stock_update=NaN(endo_nbr,gend,MAX_nsmoo);
136
  stock_innov=NaN(M_.exo_nbr,gend,MAX_ninno);
137
  if horizon
138 139
      stock_forcst_mean= NaN(endo_nbr,horizon,MAX_nforc1);
      stock_forcst_point = NaN(endo_nbr,horizon,MAX_nforc2);
140
  end
141 142
end
if nvn
143
  stock_error = NaN(length(varobs),gend,MAX_nerro);
144 145 146 147 148 149 150 151
end
if naK
    stock_filter_step_ahead =NaN(length(options_.filter_step_ahead),endo_nbr,gend+max(options_.filter_step_ahead),MAX_naK);
end
stock_param = NaN(MAX_nruns,size(myinputs.x,2));
stock_logpo = NaN(MAX_nruns,1);
stock_ys = NaN(MAX_nruns,endo_nbr);

Marco Ratto's avatar
Marco Ratto committed
152
for b=fpar:B
153
    if strcmpi(type,'prior')
154

155
        [deep, logpo] = GetOneDraw(type);
156

Marco Ratto's avatar
Marco Ratto committed
157 158
    else
        deep = x(b,:);
159 160 161 162 163
        if strcmpi(type,'posterior')
            logpo = logpost(b);
        else
            logpo = evaluate_posterior_kernel(deep');
        end
Marco Ratto's avatar
Marco Ratto committed
164
    end
165
    M_ = set_all_parameters(deep,estim_params_,M_);
166

Marco Ratto's avatar
Marco Ratto committed
167
    if run_smoother
168
        [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
Marco Ratto's avatar
Marco Ratto committed
169 170
        [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK] = ...
            DsgeSmoother(deep,gend,Y,data_index,missing_value);
171

Marco Ratto's avatar
Marco Ratto committed
172 173
        if options_.loglinear
            stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ...
174
                repmat(log(SteadyState(dr.order_var)),1,gend);
Marco Ratto's avatar
Marco Ratto committed
175
            stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ...
176
                repmat(log(SteadyState(dr.order_var)),1,gend);
Marco Ratto's avatar
Marco Ratto committed
177 178
        else
            stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ...
179
                repmat(SteadyState(dr.order_var),1,gend);
Marco Ratto's avatar
Marco Ratto committed
180
            stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ...
181
                repmat(SteadyState(dr.order_var),1,gend);
Marco Ratto's avatar
Marco Ratto committed
182 183 184 185 186 187 188 189
        end
        stock_innov(:,:,irun(2))  = etahat;
        if nvn
            stock_error(:,:,irun(3))  = epsilonhat;
        end
        if naK
            stock_filter_step_ahead(:,dr.order_var,:,irun(4)) = aK(options_.filter_step_ahead,1:endo_nbr,:);
        end
190

Marco Ratto's avatar
Marco Ratto committed
191 192 193
        if horizon
            yyyy = alphahat(iendo,i_last_obs);
            yf = forcst2a(yyyy,dr,zeros(horizon,exo_nbr));
Stéphane Adjemian's avatar
Stéphane Adjemian committed
194
            if options_.prefilter
Marco Ratto's avatar
Marco Ratto committed
195
                yf(:,IdObs) = yf(:,IdObs)+repmat(bayestopt_.mean_varobs', ...
196
                                                 horizon+maxlag,1);
Marco Ratto's avatar
Marco Ratto committed
197 198
            end
            yf(:,IdObs) = yf(:,IdObs)+(gend+[1-maxlag:horizon]')*trend_coeff';
Stéphane Adjemian's avatar
Stéphane Adjemian committed
199
            if options_.loglinear
Marco Ratto's avatar
Marco Ratto committed
200 201 202 203 204 205 206 207 208 209
                yf = yf+repmat(log(SteadyState'),horizon+maxlag,1);
            else
                yf = yf+repmat(SteadyState',horizon+maxlag,1);
            end
            yf1 = forcst2(yyyy,horizon,dr,1);
            if options_.prefilter == 1
                yf1(:,IdObs,:) = yf1(:,IdObs,:)+ ...
                    repmat(bayestopt_.mean_varobs',[horizon+maxlag,1,1]);
            end
            yf1(:,IdObs,:) = yf1(:,IdObs,:)+repmat((gend+[1-maxlag:horizon]')* ...
210
                                                   trend_coeff',[1,1,1]);
Stéphane Adjemian's avatar
Stéphane Adjemian committed
211
            if options_.loglinear
Marco Ratto's avatar
Marco Ratto committed
212 213 214 215
                yf1 = yf1 + repmat(log(SteadyState'),[horizon+maxlag,1,1]);
            else
                yf1 = yf1 + repmat(SteadyState',[horizon+maxlag,1,1]);
            end
216

217 218
            stock_forcst_mean(:,:,irun(6)) = yf(maxlag+1:end,:)';
            stock_forcst_point(:,:,irun(7)) = yf1(maxlag+1:end,:)';
Marco Ratto's avatar
Marco Ratto committed
219
        end
220 221
    else
        [T,R,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
Marco Ratto's avatar
Marco Ratto committed
222 223 224 225
    end
    stock_param(irun(5),:) = deep;
    stock_logpo(irun(5),1) = logpo;
    stock_ys(irun(5),:) = SteadyState';
226

Marco Ratto's avatar
Marco Ratto committed
227
    irun = irun +  ones(7,1);
228 229


230
    if run_smoother && (irun(1) > MAX_nsmoo || b == B),
Marco Ratto's avatar
Marco Ratto committed
231 232 233 234 235 236 237
        stock = stock_smooth(:,:,1:irun(1)-1);
        ifil(1) = ifil(1) + 1;
        save([DirectoryName '/' M_.fname '_smooth' int2str(ifil(1)) '.mat'],'stock');

        stock = stock_update(:,:,1:irun(1)-1);
        save([DirectoryName '/' M_.fname '_update' int2str(ifil(1)) '.mat'],'stock');
        if RemoteFlag==1,
238 239
            OutputFileName_smooth = [OutputFileName_smooth; {[DirectoryName filesep], [M_.fname '_smooth' int2str(ifil(1)) '.mat']}];
            OutputFileName_update = [OutputFileName_update; {[DirectoryName filesep], [M_.fname '_update' int2str(ifil(1)) '.mat']}];
Marco Ratto's avatar
Marco Ratto committed
240 241 242
        end
        irun(1) = 1;
    end
243

244
    if run_smoother && (irun(2) > MAX_ninno || b == B)
Marco Ratto's avatar
Marco Ratto committed
245 246 247 248
        stock = stock_innov(:,:,1:irun(2)-1);
        ifil(2) = ifil(2) + 1;
        save([DirectoryName '/' M_.fname '_inno' int2str(ifil(2)) '.mat'],'stock');
        if RemoteFlag==1,
249
            OutputFileName_inno = [OutputFileName_inno; {[DirectoryName filesep], [M_.fname '_inno' int2str(ifil(2)) '.mat']}];
Marco Ratto's avatar
Marco Ratto committed
250 251 252
        end
        irun(2) = 1;
    end
253

254
    if run_smoother && nvn && (irun(3) > MAX_nerro || b == B)
Marco Ratto's avatar
Marco Ratto committed
255 256 257 258
        stock = stock_error(:,:,1:irun(3)-1);
        ifil(3) = ifil(3) + 1;
        save([DirectoryName '/' M_.fname '_error' int2str(ifil(3)) '.mat'],'stock');
        if RemoteFlag==1,
259
            OutputFileName_error = [OutputFileName_error; {[DirectoryName filesep], [M_.fname '_error' int2str(ifil(3)) '.mat']}];
Marco Ratto's avatar
Marco Ratto committed
260 261 262
        end
        irun(3) = 1;
    end
263

264
    if run_smoother && naK && (irun(4) > MAX_naK || b == B)
Marco Ratto's avatar
Marco Ratto committed
265 266 267 268
        stock = stock_filter_step_ahead(:,:,:,1:irun(4)-1);
        ifil(4) = ifil(4) + 1;
        save([DirectoryName '/' M_.fname '_filter_step_ahead' int2str(ifil(4)) '.mat'],'stock');
        if RemoteFlag==1,
269
            OutputFileName_filter_step_ahead = [OutputFileName_filter_step_ahead; {[DirectoryName filesep], [M_.fname '_filter_step_ahead' int2str(ifil(4)) '.mat']}];
Marco Ratto's avatar
Marco Ratto committed
270 271 272
        end
        irun(4) = 1;
    end
273

Marco Ratto's avatar
Marco Ratto committed
274 275
    if irun(5) > MAX_nruns || b == B
        stock = stock_param(1:irun(5)-1,:);
276 277
        stock_logpo = stock_logpo(1:irun(5)-1);
        stock_ys = stock_ys(1:irun(5)-1,:);
Marco Ratto's avatar
Marco Ratto committed
278 279 280
        ifil(5) = ifil(5) + 1;
        save([DirectoryName '/' M_.fname '_param' int2str(ifil(5)) '.mat'],'stock','stock_logpo','stock_ys');
        if RemoteFlag==1,
281
            OutputFileName_param = [OutputFileName_param; {[DirectoryName filesep], [M_.fname '_param' int2str(ifil(5)) '.mat']}];
Marco Ratto's avatar
Marco Ratto committed
282 283 284
        end
        irun(5) = 1;
    end
285

286
    if run_smoother && horizon && (irun(6) > MAX_nforc1 || b == B)
Marco Ratto's avatar
Marco Ratto committed
287 288 289 290
        stock = stock_forcst_mean(:,:,1:irun(6)-1);
        ifil(6) = ifil(6) + 1;
        save([DirectoryName '/' M_.fname '_forc_mean' int2str(ifil(6)) '.mat'],'stock');
        if RemoteFlag==1,
291
            OutputFileName_forc_mean = [OutputFileName_forc_mean; {[DirectoryName filesep], [M_.fname '_forc_mean' int2str(ifil(6)) '.mat']}];
Marco Ratto's avatar
Marco Ratto committed
292 293 294
        end
        irun(6) = 1;
    end
295

296
    if run_smoother && horizon && (irun(7) > MAX_nforc2 ||  b == B)
Marco Ratto's avatar
Marco Ratto committed
297 298 299 300
        stock = stock_forcst_point(:,:,1:irun(7)-1);
        ifil(7) = ifil(7) + 1;
        save([DirectoryName '/' M_.fname '_forc_point' int2str(ifil(7)) '.mat'],'stock');
        if RemoteFlag==1,
301
            OutputFileName_forc_point = [OutputFileName_forc_point; {[DirectoryName filesep], [M_.fname '_forc_point' int2str(ifil(7)) '.mat']}];
Marco Ratto's avatar
Marco Ratto committed
302 303 304
        end
        irun(7) = 1;
    end
305
    dyn_waitbar((b-fpar+1)/(B-fpar+1),h);
Marco Ratto's avatar
Marco Ratto committed
306 307 308 309
end

myoutput.ifil=ifil;
if RemoteFlag==1,
310
    myoutput.OutputFileName = [OutputFileName_smooth;
311 312 313 314 315 316 317
                        OutputFileName_update;
                        OutputFileName_inno;
                        OutputFileName_error;
                        OutputFileName_filter_step_ahead;
                        OutputFileName_param;
                        OutputFileName_forc_mean;
                        OutputFileName_forc_point];
Marco Ratto's avatar
Marco Ratto committed
318
end
319

320 321
dyn_waitbar_close(h);