stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fixes with screen output for parallel/octave + cosmethics

time-shift
Marco Ratto 2011-02-02 14:13:11 +01:00
parent 8111b671ff
commit bed32115d2
5 changed files with 171 additions and 147 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
matlab/McMCDiagnostics_core.m
View File

@ -2,17 +2,17 @@ function myoutput = McMCDiagnostics_core(myinputs,fpar,npar,whoiam, ThisMatlab)
% PARALLEL CONTEXT
% Core functionality for MCMC Diagnostics, which can be parallelized.
% See also the comment in random_walk_metropolis_hastings_core.m funtion.
% INPUTS
% INPUTS
% See See the comment in random_walk_metropolis_hastings_core.m funtion.
% OUTPUTS
% o myoutput [struc]
% Contained UDIAG.
%
% ALGORITHM
% Portion of McMCDiagnostics.m function.
% ALGORITHM
% Portion of McMCDiagnostics.m function.
%
% SPECIAL REQUIREMENTS.
% None.
@ -59,7 +59,7 @@ if ~exist('MhDirectoryName'),
MhDirectoryName = CheckPath('metropolis');
end
ALPHA = 0.2; % increase too much with the number of simulations.
ALPHA = 0.2; % increase too much with the number of simulations.
tmp = zeros(NumberOfDraws*nblck,3);
UDIAG = zeros(NumberOfLines,6,npar-fpar+1);
@ -69,30 +69,36 @@ if whoiam
waitbarTitle=['Local '];
else
waitbarTitle=[Parallel(ThisMatlab).ComputerName];
end
fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab));
end
fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab));
end
for j=fpar:npar,
fprintf(' Parameter %d... ',j);
if exist('OCTAVE_VERSION'),
if (whoiam==0),
printf(' Parameter %d... ',j);
end
else
fprintf(' Parameter %d... ',j);
end
for b = 1:nblck
startline = 0;
for n = 1:NumberOfMcFilesPerBlock
%load([MhDirectoryName '/' mcfiles(n,1,b).name],'x2');
load([MhDirectoryName '/' M_.fname '_mh',int2str(n),'_blck' int2str(b) ...
'.mat'],'x2');
'.mat'],'x2');
nx2 = size(x2,1);
tmp((b-1)*NumberOfDraws+startline+(1:nx2),1) = x2(:,j);
% clear x2;
startline = startline + nx2;
end
% $$$ %load([MhDirectoryName '/' mcfiles(NumberOfMcFilesPerBlock,1,b).name],'x2');
% $$$ load([MhDirectoryName '/' M_.fname '_mh',int2str(NumberOfMcFilesPerBlock),'_blck' int2str(b) '.mat'],'x2');
% $$$ tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+MAX_nruns*(LastFileNumber-1)+LastLineNumber,1) = x2(:,j);
% $$$ clear x2;
% $$$ startline = startline + LastLineNumber;
% $$$ %load([MhDirectoryName '/' mcfiles(NumberOfMcFilesPerBlock,1,b).name],'x2');
% $$$ load([MhDirectoryName '/' M_.fname '_mh',int2str(NumberOfMcFilesPerBlock),'_blck' int2str(b) '.mat'],'x2');
% $$$ tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+MAX_nruns*(LastFileNumber-1)+LastLineNumber,1) = x2(:,j);
% $$$ clear x2;
% $$$ startline = startline + LastLineNumber;
end
tmp(:,2) = kron(transpose(1:nblck),ones(NumberOfDraws,1));
tmp(:,3) = kron(ones(nblck,1),time');
tmp(:,3) = kron(ones(nblck,1),time');
tmp = sortrows(tmp,1);
ligne = 0;
for iter = Origin:StepSize:NumberOfDraws
@ -111,13 +117,19 @@ for j=fpar:npar,
for i=1:nblck
pmet = temp(find(temp(:,2)==i));
UDIAG(ligne,2,j-fpar+1) = UDIAG(ligne,2,j-fpar+1) + pmet(csup,1)-pmet(cinf,1);
moyenne = mean(pmet,1); %% Within mean.
moyenne = mean(pmet,1); %% Within mean.
UDIAG(ligne,4,j-fpar+1) = UDIAG(ligne,4,j-fpar+1) + sum((pmet(:,1)-moyenne).^2)/(n-1);
UDIAG(ligne,6,j-fpar+1) = UDIAG(ligne,6,j-fpar+1) + sum(abs(pmet(:,1)-moyenne).^3)/(n-1);
end
end
fprintf('Done! \n');
if whoiam,
if exist('OCTAVE_VERSION'),
if (whoiam==0),
printf('Done! \n');
end
else
fprintf('Done! \n');
end
if whoiam,
waitbarString = [ 'Parameter ' int2str(j) '/' int2str(npar) ' done.'];
fMessageStatus((j-fpar+1)/(npar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab))
end

72
matlab/PosteriorIRF_core1.m
View File

@ -5,10 +5,10 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,npar,whoiam, ThisMatlab)
% that running in parallel a 'for' cycle, this function run in parallel a
% 'while' loop! The parallelization of 'while' loop (when possible) is a more
% sophisticated procedure.
%
%
% See also the comment in random_walk_metropolis_hastings_core.m funtion.
%
% INPUTS
% INPUTS
% See the comment in random_walk_metropolis_hastings_core.m funtion.
%
% OUTPUTS
@ -16,8 +16,8 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,npar,whoiam, ThisMatlab)
% Contained:
% OutputFileName_dsge, OutputFileName_param and OutputFileName_bvardsge.
%
% ALGORITHM
% Portion of PosteriorIRF.m function. Specifically the 'while' cycle.
% ALGORITHM
% Portion of PosteriorIRF.m function. Specifically the 'while' cycle.
%
% SPECIAL REQUIREMENTS.
% None.
@ -109,10 +109,10 @@ else
h = waitbar(0,'Bayesian (prior) IRFs...');
end
end
end
OutputFileName_bvardsge = {};
OutputFileName_bvardsge = {};
OutputFileName_dsge = {};
OutputFileName_param = {};
@ -121,25 +121,25 @@ fpar0=fpar;
fpar = fpar-1;
if whoiam
ifil2=ifil2(whoiam);
NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge(whoiam);
NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar(whoiam);
ifil2=ifil2(whoiam);
NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge(whoiam);
NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar(whoiam);
end
% Parallel 'while' very good!!!
while fpar<npar
while fpar<npar
fpar = fpar + 1;
irun = irun+1;
irun2 = irun2+1;
if strcmpi(type,'prior')
if strcmpi(type,'prior')
deep = GetOneDraw(type);
else
deep = x(fpar,:);
end
stock_param(irun2,:) = deep;
stock_param(irun2,:) = deep;
set_parameters(deep);
[dr,info] = resol(oo_.steady_state,0);
[dr,info] = resol(oo_.steady_state,0);
if info(1)
nosaddle = nosaddle + 1;
fpar = fpar - 1;
@ -168,7 +168,7 @@ while fpar<npar
y = 100*y/cs(i,i);
end
for j = 1:nvar
if max(y(IndxVariables(j),:)) - min(y(IndxVariables(j),:)) > 1e-12
if max(y(IndxVariables(j),:)) - min(y(IndxVariables(j),:)) > 1e-12
stock_irf_dsge(:,j,i,irun) = transpose(y(IndxVariables(j),:));
end
end
@ -179,21 +179,21 @@ while fpar<npar
[fval,cost_flag,info,PHI,SIGMAu,iXX] = DsgeVarLikelihood(deep',gend);
dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names));
DSGE_PRIOR_WEIGHT = floor(gend*(1+dsge_prior_weight));
SIGMA_inv_upper_chol = chol(inv(SIGMAu*gend*(dsge_prior_weight+1)));
SIGMA_inv_upper_chol = chol(inv(SIGMAu*gend*(dsge_prior_weight+1)));
explosive_var = 1;
while explosive_var
% draw from the marginal posterior of SIGMA
SIGMAu_draw = rand_inverse_wishart(nvobs, DSGE_PRIOR_WEIGHT-NumberOfParametersPerEquation, ...
SIGMA_inv_upper_chol);
SIGMA_inv_upper_chol);
% draw from the conditional posterior of PHI
PHI_draw = rand_matrix_normal(NumberOfParametersPerEquation,nvobs, PHI, ...
chol(SIGMAu_draw)', chol(iXX)');
chol(SIGMAu_draw)', chol(iXX)');
Companion_matrix(1:nvobs,:) = transpose(PHI_draw(1:NumberOfLagsTimesNvobs,:));
% Check for stationarity
explosive_var = any(abs(eig(Companion_matrix))>1.000000001);
end
% Get the mean
mu = zeros(1,nvobs);
% Get the mean
mu = zeros(1,nvobs);
% Get rotation
if dsge_prior_weight > 0
Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e);
@ -213,7 +213,7 @@ while fpar<npar
end
tmp_dsgevar = kron(ones(options_.irf,1),mu);
for j = 1:(nvobs*M_.exo_nbr)
if max(irfs(:,j)) - min(irfs(:,j)) > 1e-10
if max(irfs(:,j)) - min(irfs(:,j)) > 1e-10
tmp_dsgevar(:,j) = (irfs(:,j));
end
end
@ -222,12 +222,12 @@ while fpar<npar
else
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,nvobs,M_.exo_nbr);
instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' ...
int2str(NumberOfIRFfiles_dsgevar) '.mat stock_irf_bvardsge;'];,
int2str(NumberOfIRFfiles_dsgevar) '.mat stock_irf_bvardsge;'];,
eval(['save ' instr]);
if RemoteFlag==1,
OutputFileName_bvardsge = [OutputFileName_bvardsge; {[MhDirectoryName filesep], [M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) '.mat']}];
end
NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
IRUN =0;
stock_irf_dsgevar = zeros(options_.irf,nvobs,M_.exo_nbr,MAX_nirfs_dsgevar);
end
@ -238,7 +238,7 @@ while fpar<npar
if MAX_nirfs_dsgevar & (fpar == npar | IRUN == npar)
stock_irf_bvardsge = stock_irf_bvardsge(:,:,:,1:IRUN);
instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' ...
int2str(NumberOfIRFfiles_dsgevar) '.mat stock_irf_bvardsge;'];,
int2str(NumberOfIRFfiles_dsgevar) '.mat stock_irf_bvardsge;'];,
eval(['save ' instr]);
NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
if RemoteFlag==1,
@ -250,7 +250,7 @@ while fpar<npar
save([MhDirectoryName '/' M_.fname '_irf_dsge' int2str(NumberOfIRFfiles_dsge) '.mat'],'stock_irf_dsge');
if RemoteFlag==1,
OutputFileName_dsge = [OutputFileName_dsge; {[MhDirectoryName filesep], [M_.fname '_irf_dsge' int2str(NumberOfIRFfiles_dsge) '.mat']}];
end
end
NumberOfIRFfiles_dsge = NumberOfIRFfiles_dsge+1;
irun = 0;
end
@ -266,21 +266,25 @@ while fpar<npar
ifil2 = ifil2 + 1;
irun2 = 0;
end
if exist('OCTAVE_VERSION')
printf(['Posterior IRF %3.f%% done\r'],(fpar/npar*100));
elseif ~whoiam
waitbar(fpar/npar,h);
if exist('OCTAVE_VERSION'),
if (whoiam==0),
printf(['Posterior IRF %3.f%% done\r'],(fpar/npar*100));
end
elseif ~whoiam,
waitbar(fpar/npar,h);
end
if whoiam,
if ~exist('OCTAVE_VERSION')
fprintf('Done! \n');
end
waitbarString = [ 'Subdraw ' int2str(fpar) '/' int2str(npar) ' done.'];
fMessageStatus((fpar-fpar0)/(npar-fpar0),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab));
end
if whoiam,
fprintf('Done! \n');
waitbarString = [ 'Subdraw ' int2str(fpar) '/' int2str(npar) ' done.'];
fMessageStatus((fpar-fpar0)/(npar-fpar0),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab));
end
end
if whoiam==0
if nosaddle
disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(npar+nosaddle))])
disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(npar+nosaddle))])
end
if exist('h')
close(h);

32
matlab/independent_metropolis_hastings_core.m
View File

@ -3,15 +3,15 @@ function myoutput = independent_metropolis_hastings_core(myinputs,fblck,nblck,wh
% The most computationally intensive portion of code in
% independent_metropolis_hastings (the 'for xxx = fblck:nblck' cycle).
% See the comment in random_walk_metropolis_hastings_core.m funtion.
%
% INPUTS
%
% INPUTS
% See See the comment in random_walk_metropolis_hastings_core.m funtion.
% OUTPUTS
% See See the comment in random_walk_metropolis_hastings_core.m funtion.
%
% ALGORITHM
% Portion of Independing Metropolis-Hastings.
% ALGORITHM
% Portion of Independing Metropolis-Hastings.
%
% SPECIAL REQUIREMENTS.
% None.
@ -105,7 +105,7 @@ for b = fblck:nblck,
rand('state',record.Seeds(b).Unifor);
if (options_.load_mh_file~=0) & (fline(b)>1) & OpenOldFile(b)
load(['./' MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) ...
'_blck' int2str(b) '.mat'])
'_blck' int2str(b) '.mat'])
x2 = [x2;zeros(InitSizeArray(b)-fline(b)+1,npar)];
logpo2 = [logpo2;zeros(InitSizeArray(b)-fline(b)+1,1)];
OpenOldFile(b) = 0;
@ -124,22 +124,22 @@ for b = fblck:nblck,
waitbarTitle=['Local '];
else
waitbarTitle=[options_.parallel(ThisMatlab).ComputerName];
end
end
fMessageStatus(0,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab));
else,
hh = waitbar(0,['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...']);
set(hh,'Name','Metropolis-Hastings');
end
isux = 0;
jsux = 0;
irun = fline(b);
j = 1;
while j <= nruns(b)
par = feval(ProposalFun, xparam1, proposal_covariance, n);
par = feval(ProposalFun, xparam1, proposal_covariance, n);
if all( par(:) > mh_bounds(:,1) ) & all( par(:) < mh_bounds(:,2) )
try
logpost = - feval(TargetFun, par(:),varargin{:});
logpost = - feval(TargetFun, par(:),varargin{:});
catch,
logpost = -inf;
end
@ -152,11 +152,11 @@ for b = fblck:nblck,
if (logpost > -inf) && (log(rand) < r)
x2(irun,:) = par;
ix2(b,:) = par;
logpo2(irun) = logpost;
logpo2(irun) = logpost;
ilogpo2(b) = logpost;
isux = isux + 1;
jsux = jsux + 1;
else
else
x2(irun,:) = ix2(b,:);
logpo2(irun) = ilogpo2(b);
end
@ -164,12 +164,14 @@ for b = fblck:nblck,
if exist('OCTAVE_VERSION') || options_.console_mode
if mod(j, 10) == 0
if exist('OCTAVE_VERSION')
printf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acception rate: %3.f%%\r', b, nblck, 100 * prtfrc, 100 * isux / j);
if (whoiam==0),
printf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acception rate: %3.f%%\r', b, nblck, 100 * prtfrc, 100 * isux / j);
end
else
fprintf(' MH: Computing Metropolis-Hastings (chain %d/%d): %3.f \b%% done, acception rate: %3.f \b%%\r', b, nblck, 100 * prtfrc, 100 * isux / j);
end
end
if mod(j,50)==0 & whoiam,
if mod(j,50)==0 & whoiam,
% keyboard;
waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) '), ' sprintf('accept. %3.f%%%%', 100 * isux/j)];
fMessageStatus(prtfrc,whoiam,waitbarString, '', options_.parallel(ThisMatlab))
@ -177,13 +179,13 @@ for b = fblck:nblck,
else
if mod(j, 3)==0 & ~whoiam
waitbar(prtfrc,hh,[ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]);
elseif mod(j,50)==0 & whoiam,
elseif mod(j,50)==0 & whoiam,
% keyboard;
waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)];
fMessageStatus(prtfrc,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab))
end
end
if (irun == InitSizeArray(b)) | (j == nruns(b)) % Now I save the simulations
save([MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'],'x2','logpo2');
fidlog = fopen([MhDirectoryName '/metropolis.log'],'a');

126
matlab/prior_posterior_statistics_core.m
View File

@ -3,22 +3,22 @@ function myoutput=prior_posterior_statistics_core(myinputs,fpar,B,whoiam, ThisMa
% Core functionality for prior_posterior.m function, which can be parallelized.
% See also the comment in random_walk_metropolis_hastings_core.m funtion.
%
% INPUTS
% 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;
% Contained OutputFileName_smooth;
% _update;
% _inno;
% _error;
% _filter_step_ahead;
% _param;
% _forc_mean;
% _forc_point
%
% ALGORITHM
% Portion of prior_posterior.m function.
% ALGORITHM
% Portion of prior_posterior.m function.
% This file is part of Dynare.
%
% SPECIAL REQUIREMENTS.
@ -111,37 +111,37 @@ else
else
h = waitbar(0,'Taking subdraws...');
end
end
if RemoteFlag==1,
OutputFileName_smooth = {};
OutputFileName_update = {};
OutputFileName_inno = {};
OutputFileName_error = {};
OutputFileName_filter_step_ahead = {};
OutputFileName_param = {};
OutputFileName_forc_mean = {};
OutputFileName_forc_point = {};
% OutputFileName_moments = {};
OutputFileName_smooth = {};
OutputFileName_update = {};
OutputFileName_inno = {};
OutputFileName_error = {};
OutputFileName_filter_step_ahead = {};
OutputFileName_param = {};
OutputFileName_forc_mean = {};
OutputFileName_forc_point = {};
% OutputFileName_moments = {};
end
for b=fpar:B
% [deep, logpo] = GetOneDraw(type);
% set_all_parameters(deep);
% dr = resol(oo_.steady_state,0);
if strcmpi(type,'prior')
[deep, logpo] = GetOneDraw(type);
else
deep = x(b,:);
logpo = logpost(b);
end
set_all_parameters(deep);
[dr,info] = resol(oo_.steady_state,0);
if run_smoother
[alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK] = ...
DsgeSmoother(deep,gend,Y,data_index,missing_value);
@ -164,7 +164,7 @@ for b=fpar:B
if naK
stock_filter_step_ahead(:,dr.order_var,:,irun(4)) = aK(options_.filter_step_ahead,1:endo_nbr,:);
end
if horizon
yyyy = alphahat(iendo,i_last_obs);
yf = forcst2a(yyyy,dr,zeros(horizon,exo_nbr));
@ -190,19 +190,19 @@ for b=fpar:B
else
yf1 = yf1 + repmat(SteadyState',[horizon+maxlag,1,1]);
end
stock_forcst_mean(:,:,irun(6)) = yf';
stock_forcst_point(:,:,irun(7)) = yf1';
end
end
stock_param(irun(5),:) = deep;
stock_logpo(irun(5),1) = logpo;
stock_ys(irun(5),:) = SteadyState';
irun = irun + ones(7,1);
if irun(1) > MAX_nsmoo || b == B
stock = stock_smooth(:,:,1:irun(1)-1);
ifil(1) = ifil(1) + 1;
@ -211,72 +211,72 @@ for b=fpar:B
stock = stock_update(:,:,1:irun(1)-1);
save([DirectoryName '/' M_.fname '_update' int2str(ifil(1)) '.mat'],'stock');
if RemoteFlag==1,
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']}];
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']}];
end
irun(1) = 1;
end
if irun(2) > MAX_ninno || b == B
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,
OutputFileName_inno = [OutputFileName_inno; {[DirectoryName filesep], [M_.fname '_inno' int2str(ifil(2)) '.mat']}];
OutputFileName_inno = [OutputFileName_inno; {[DirectoryName filesep], [M_.fname '_inno' int2str(ifil(2)) '.mat']}];
end
irun(2) = 1;
end
if nvn && (irun(3) > MAX_nerro || b == B)
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,
OutputFileName_error = [OutputFileName_error; {[DirectoryName filesep], [M_.fname '_error' int2str(ifil(3)) '.mat']}];
OutputFileName_error = [OutputFileName_error; {[DirectoryName filesep], [M_.fname '_error' int2str(ifil(3)) '.mat']}];
end
irun(3) = 1;
end
if naK && (irun(4) > MAX_naK || b == B)
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,
OutputFileName_filter_step_ahead = [OutputFileName_filter_step_ahead; {[DirectoryName filesep], [M_.fname '_filter_step_ahead' int2str(ifil(4)) '.mat']}];
OutputFileName_filter_step_ahead = [OutputFileName_filter_step_ahead; {[DirectoryName filesep], [M_.fname '_filter_step_ahead' int2str(ifil(4)) '.mat']}];
end
irun(4) = 1;
end
if irun(5) > MAX_nruns || b == B
stock = stock_param(1:irun(5)-1,:);
ifil(5) = ifil(5) + 1;
save([DirectoryName '/' M_.fname '_param' int2str(ifil(5)) '.mat'],'stock','stock_logpo','stock_ys');
if RemoteFlag==1,
OutputFileName_param = [OutputFileName_param; {[DirectoryName filesep], [M_.fname '_param' int2str(ifil(5)) '.mat']}];
OutputFileName_param = [OutputFileName_param; {[DirectoryName filesep], [M_.fname '_param' int2str(ifil(5)) '.mat']}];
end
irun(5) = 1;
end
if horizon && (irun(6) > MAX_nforc1 || b == B)
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,
OutputFileName_forc_mean = [OutputFileName_forc_mean; {[DirectoryName filesep], [M_.fname '_forc_mean' int2str(ifil(6)) '.mat']}];
OutputFileName_forc_mean = [OutputFileName_forc_mean; {[DirectoryName filesep], [M_.fname '_forc_mean' int2str(ifil(6)) '.mat']}];
end
irun(6) = 1;
end
if horizon && (irun(7) > MAX_nforc2 || b == B)
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,
OutputFileName_forc_point = [OutputFileName_forc_point; {[DirectoryName filesep], [M_.fname '_forc_point' int2str(ifil(7)) '.mat']}];
OutputFileName_forc_point = [OutputFileName_forc_point; {[DirectoryName filesep], [M_.fname '_forc_point' int2str(ifil(7)) '.mat']}];
end
irun(7) = 1;
end
% if moments_varendo && (irun(8) > MAX_momentsno || b == B)
% stock = stock_moments(1:irun(8)-1);
% ifil(8) = ifil(8) + 1;
@ -286,18 +286,22 @@ for b=fpar:B
% end
% irun(8) = 1;
% end
% DirectoryName=TempPath;
if exist('OCTAVE_VERSION')
printf('Taking subdraws: %3.f%% done\r', b/B*100);
% DirectoryName=TempPath;
if exist('OCTAVE_VERSION'),
if (whoiam==0),
printf('Taking subdraws: %3.f%% done\r', b/B*100);
end
elseif ~whoiam,
waitbar(b/B,h);
end
if whoiam,
fprintf('Done! \n');
if ~exist('OCTAVE_VERSION')
fprintf('Done! \n');
end
waitbarString = [ 'Subdraw ' int2str(b) '/' int2str(B) ' done.'];
fMessageStatus((b-fpar+1)/(B-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab));
end
@ -305,15 +309,15 @@ end
myoutput.ifil=ifil;
if RemoteFlag==1,
myoutput.OutputFileName = [OutputFileName_smooth;
OutputFileName_update;
OutputFileName_inno;
OutputFileName_error;
OutputFileName_filter_step_ahead;
OutputFileName_param;
OutputFileName_forc_mean;
OutputFileName_forc_point];
% OutputFileName_moments];
myoutput.OutputFileName = [OutputFileName_smooth;
OutputFileName_update;
OutputFileName_inno;
OutputFileName_error;
OutputFileName_filter_step_ahead;
OutputFileName_param;
OutputFileName_forc_mean;
OutputFileName_forc_point];
% OutputFileName_moments];
end
if exist('OCTAVE_VERSION')
@ -323,5 +327,5 @@ else
if exist('h')
close(h)
end
end

40
matlab/random_walk_metropolis_hastings_core.m
View File

@ -3,8 +3,8 @@ function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,wh
% This function contain the most computationally intensive portion of code in
% random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop). The branches in 'for'
% cycle and are completely independent than suitable to be executed in parallel way.
%
% INPUTS
%
% INPUTS
% o myimput [struc] The mandatory variables for local/remote
% parallel computing obtained from random_walk_metropolis_hastings.m
% function.
@ -26,8 +26,8 @@ function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,wh
% NewFile;
% OutputFileName
%
% ALGORITHM
% Portion of Metropolis-Hastings.
% ALGORITHM
% Portion of Metropolis-Hastings.
%
% SPECIAL REQUIREMENTS.
% None.
@ -88,16 +88,16 @@ nruns=myinputs.nruns;
NewFile=myinputs.NewFile;
MAX_nruns=myinputs.MAX_nruns;
d=myinputs.d;
InitSizeArray=myinputs.InitSizeArray;
InitSizeArray=myinputs.InitSizeArray;
record=myinputs.record;
varargin=myinputs.varargin;
% Necessary only for remote computing!
if whoiam
Parallel=myinputs.Parallel;
% initialize persistent variables in priordens()
priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, ...
bayestopt_.p3,bayestopt_.p4,1);
Parallel=myinputs.Parallel;
% initialize persistent variables in priordens()
priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, ...
bayestopt_.p3,bayestopt_.p4,1);
end
% (re)Set the penalty
@ -138,7 +138,7 @@ for b = fblck:nblck,
rand('state',record.Seeds(b).Unifor);
if (options_.load_mh_file~=0) & (fline(b)>1) & OpenOldFile(b)
load(['./' MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) ...
'_blck' int2str(b) '.mat'])
'_blck' int2str(b) '.mat'])
x2 = [x2;zeros(InitSizeArray(b)-fline(b)+1,npar)];
logpo2 = [logpo2;zeros(InitSizeArray(b)-fline(b)+1,1)];
OpenOldFile(b) = 0;
@ -160,14 +160,14 @@ for b = fblck:nblck,
hh = waitbar(0,['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...']);
set(hh,'Name','Metropolis-Hastings');
end
end
if whoiam
if options_.parallel(ThisMatlab).Local,
waitbarTitle=['Local '];
else
waitbarTitle=[options_.parallel(ThisMatlab).ComputerName];
end
end
prc0=(b-fblck)/(nblck-fblck+1)*(exist('OCTAVE_VERSION') || options_.console_mode);
fMessageStatus(prc0,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab));
end
@ -179,7 +179,7 @@ for b = fblck:nblck,
par = feval(ProposalFun, ix2(b,:), proposal_covariance_Cholesky_decomposition, n);
if all( par(:) > mh_bounds(:,1) ) & all( par(:) < mh_bounds(:,2) )
try
logpost = - feval(TargetFun, par(:),varargin{:});
logpost = - feval(TargetFun, par(:),varargin{:});
catch
logpost = -inf;
end
@ -189,11 +189,11 @@ for b = fblck:nblck,
if (logpost > -inf) && (log(rand) < logpost-ilogpo2(b))
x2(irun,:) = par;
ix2(b,:) = par;
logpo2(irun) = logpost;
logpo2(irun) = logpost;
ilogpo2(b) = logpost;
isux = isux + 1;
jsux = jsux + 1;
else
else
x2(irun,:) = ix2(b,:);
logpo2(irun) = ilogpo2(b);
end
@ -201,14 +201,16 @@ for b = fblck:nblck,
if exist('OCTAVE_VERSION') || options_.console_mode
if mod(j, 10) == 0
if exist('OCTAVE_VERSION')
printf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acception rate: %3.f%%\r', b, nblck, 100 * prtfrc, 100 * isux / j);
if (whoiam==0)
printf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acception rate: %3.f%%\r', b, nblck, 100 * prtfrc, 100 * isux / j);
end
else
s0=repmat('\b',1,length(newString));
newString=sprintf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acceptance rate: %3.f%%', b, nblck, 100 * prtfrc, 100 * isux / j);
fprintf([s0,'%s'],newString);
end
end
if mod(j,50)==0 & whoiam
if mod(j,50)==0 & whoiam
% keyboard;
waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) '), ' sprintf('accept. %3.f%%', 100 * isux/j)];
fMessageStatus((b-fblck)/(nblck-fblck+1)+prtfrc/(nblck-fblck+1),whoiam,waitbarString, '', options_.parallel(ThisMatlab));
@ -216,13 +218,13 @@ for b = fblck:nblck,
else
if mod(j, 3)==0 & ~whoiam
waitbar(prtfrc,hh,[ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]);
elseif mod(j,50)==0 & whoiam,
elseif mod(j,50)==0 & whoiam,
% keyboard;
waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)];
fMessageStatus(prtfrc,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab));
end
end
if (irun == InitSizeArray(b)) | (j == nruns(b)) % Now I save the simulations
save([MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'],'x2','logpo2');
fidlog = fopen([MhDirectoryName '/metropolis.log'],'a');