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

Merge branch 'master' of ssh://kirikou.dynare.org/srv/d_kirikou/git/dynare-ratto

time-shift
Stéphane Adjemian (Karaba) 2010-02-13 12:45:40 +01:00
parent 8fa10d61cc be510a4b93
commit 498c10b2ab
20 changed files with 1407 additions and 557 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

131
matlab/AnalyseComputationalEnviroment.m
View File

@ -1,5 +1,12 @@
function [ErrorCode] = AnalyseComputationalEnviroment(DataInput)
% Input/Output description:
% DESCRIPTION
% This function is used to check the user computational request.
% If no error happen the function return 0.
% INPUT/OUTPUT description:
%
% DataInput is the strcture option_.parallel, with the follow fields:
%
@ -37,10 +44,17 @@ function [ErrorCode] = AnalyseComputationalEnviroment(DataInput)
%
% Value 3: The remote computer is unreachable!!!
%
% Value 4: The user name and/or password is/are incorrect on the
% remote computer!
% Value 4: The fields user name and/or password are/is empty!
%
% Value 5: Remote Drive and/or Remote Folder not exist!
%
% Value 6: It is impossible write/read file on remote computer.
%
% Value 7: The valu user and/or passwd are incorret or the user have
% no permission to execute a Matlab section.
%
%
%
% Value 5: It is impossible write/read file on remote computer.
%
% Then at the point call of this function it is possible react in a best way, in accord
% with the ErrorCode.
@ -95,13 +109,13 @@ if (DataInput.Local == 1)
% We look for the information on local computer hardware.
si=[];
de=[];
si0=[];
de0=[];
[si de]=system(['psinfo \\']);
[si0 de0]=system(['psinfo \\']);
RealNumCPU=-1;
RealNumCPU=GiveCPUnumber(de);
RealNumCPU=GiveCPUnumber(de0);
% Trasforming the input data provided in a form [n1:n2] in a single numerical
% value.
@ -125,15 +139,15 @@ end
% In this case we need more sophisticated check.
if (DataInput.Local == 0)
si=[];
de=[];
si1=[];
de1=[];
[si de]=system(['ping ', DataInput.PcName]);
[si1 de1]=system(['ping ', DataInput.PcName]);
if si==1
if si1==1
% It is impossiblie to be connected to the
% remote computer.
@ -141,56 +155,93 @@ if (DataInput.Local == 0)
return;
end
% -> IL CODICE SEGUENTE E' DA CONTROLLARE E VERIFICARE!
% The Local Machine can be connetted with Remote Computer.
% Now we verify if user name and password are correct and if remote
% drive and remote folder exist on the remote computer and it is
% possible to exchange data with them.
if (isempty(DataInput.user)) || (isempty(DataInput.passwd))
% The fields user name and/or password are/is empty!
ErrorCode=4;
return
si=[];
de=[];
[si de]=system(['psinfo \\', DataInput.PcName, ' -u ',DataInput.user, ' -p ',DataInput.passwd ]);
if si<0
% It is possible to be connected to the remote computer but it is not usable because the user
% name and/or password is/are incorrect.
ErrorCodeComputer=4;
return;
else
% Username and Password are correct!
end
% Now we verify if it possible to exchange data with the remote
% computer:
% Now we very if RemoteDrive and/or RemoteFolder exist on remote
% computer
StartPwd=pwd;
fid = fopen('Tracing.txt','w+');
fclose (fid);
try
cd(['\\',DataInput.PcName,'\',DataInput.RemoteDrive,'$\',DataInput.RemoteFolder]);
catch
cd ([StartPwd]);
disp 'Remote Drive and/or Remote Folder not exist!';
ErrorCode=5;
return
end
% ATTENZIONE: verificare perche sembra funzionare anche se il RemoteFolder non
% esiste.
cd ([StartPwd]);
Status=movefile('Tracing.txt', ['\\',DataInput.PcName,'\',DataInput.RemoteDrive,'$\',DataInput.RemoteFolder]);
% Now we verify if it possible to exchange data with the remote
% computer:
Status=copyfile('Tracing.m', ['\\',DataInput.PcName,'\',DataInput.RemoteDrive,'$\',DataInput.RemoteFolder]);
if Status==1
% Remote Drive/Folder exist on Remote computer and
% it is possible to exchange data with him.
else
% Move file error!
ErrorCodeComputer=5;
ErrorCode=6;
return;
end
% Now we verify if it is possible execute a matlab section on remote
% machine when the user is .user with password .passwd
si2=[];
de2=[];
[si2 de2]=system(['start /B /WAIT psexec \\',DataInput.PcName,' -e -u ',DataInput.user,' -p ',DataInput.passwd,' -W ',DataInput.RemoteDrive,':\',DataInput.RemoteFolder, ' -low matlab -nosplash -nodesktop -minimize -r Tracing']);
NoError='error code 0';
StrError= findstr(NoError,de2);
if isempty (StrError)
% Bad user and/or passwd!
ErrorCode=7;
return;
else
% No error it is possible execute a matlab section on remote
% machine when the user is .user with password .passwd
end
% At this point we can to analyze the remote computer hardware.
RealNumCPU=-1;
RealNumCPU=GiveCPUnumber(de);
[si0 de0]=system(['psinfo \\']);
RealNumCPU=GiveCPUnumber(de0);
% Trasforming the input data provided in a form [n1:n2] in a single numerical
% value.

4
matlab/DiffuseKalmanSmoother1.m
View File

@ -131,7 +131,7 @@ while notsteady & t<smpl
a(:,t+1) = T*atilde(:,t);
aK(1,:,t+1) = a(:,t+1);
for jnk=2:nk,
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1))';
end
P(:,:,t+1) = T*P(:,:,t)*transpose(T)-T*P(:,mf,t)*transpose(K(:,:,t)) + QQ;
notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
@ -154,7 +154,7 @@ while t<smpl
a(:,t+1) = T*a(:,t) + K_s*v(:,t);
aK(1,:,t+1) = a(:,t+1);
for jnk=2:nk,
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1))';
end
end
t = smpl+1;

58
matlab/GiveCPUnumber.m Normal file
View File

@ -0,0 +1,58 @@
function [nCPU]= GiveCPUnumber (ComputerInformations)
% DESCRIPTION
% This function return the CPUs or cores numer avaiable
% on the computer used for run parallel code.
% INPUTS
% an array contained several fields that describe the hardaware
% software enviroments of a generic computer.
%
% OUTPUTS
% The CPUs or Cores numbers of computer.
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2005-2009 Dynare Team
%
% This file is part of Dynare.
%
% 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/>.
nCPU=-1;
OffSet=27;
SringPosition=strfind(ComputerInformations, 'Processors:');
nCPU=ComputerInformations(SringPosition+OffSet);
% We check if there are Processors/Cores more than 9.
t0=ComputerInformations(SringPosition+OffSet+1);
t1=str2num(t0);
t1=isempty(t1);
% if t1 is 0 the machine have more than 9 CPU.
if t1==0
nCPU=strcat(nCPU,t0);
end
nCPU=str2num(nCPU);
return

2
matlab/McMCDiagnostics.m
View File

@ -101,7 +101,7 @@ if isnumeric(options_.parallel),
UDIAG = fout.UDIAG;
clear fout
else
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, 1, npar,{},'McMCDiagnostics_core', localVars);
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, 1, npar,{},'McMCDiagnostics_core', localVars, [], options_.parallel_info);
UDIAG = fout(1).UDIAG;
for j=2:totCPU,
UDIAG = cat(3,UDIAG ,fout(j).UDIAG);

518
matlab/PosteriorIRF.m
View File

@ -28,6 +28,7 @@ function PosteriorIRF(type)
% 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_ estim_params_ oo_ M_ bayestopt_
% Set the number of periods
if isempty(options_.irf) | ~options_.irf
@ -71,32 +72,32 @@ else
end
DirectoryName = CheckPath('Output');
if strcmpi(type,'posterior')
MhDirectoryName = CheckPath('metropolis');
MhDirectoryName = CheckPath('metropolis');
elseif strcmpi(type,'gsa')
MhDirectoryName = CheckPath('GSA');
MhDirectoryName = CheckPath('GSA');
else
MhDirectoryName = CheckPath('prior');
MhDirectoryName = CheckPath('prior');
end
if strcmpi(type,'posterior')
load([ MhDirectoryName '/' M_.fname '_mh_history.mat'])
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
load([ MhDirectoryName filesep M_.fname '_mh_history.mat'])
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
elseif strcmpi(type,'gsa')
load([ MhDirectoryName '/' M_.fname '_prior.mat'],'lpmat0','lpmat','istable')
x=[lpmat0(istable,:) lpmat(istable,:)];
clear lpmat istable
NumberOfDraws=size(x,1);
B=NumberOfDraws; options_.B = B;
load([ MhDirectoryName filesep M_.fname '_prior.mat'],'lpmat0','lpmat','istable')
x=[lpmat0(istable,:) lpmat(istable,:)];
clear lpmat istable
NumberOfDraws=size(x,1);
B=NumberOfDraws; options_.B = B;
else% type = 'prior'
NumberOfDraws = 500;
NumberOfDraws = 500;
end
if ~strcmpi(type,'gsa')
B = min([round(.5*NumberOfDraws),500]); options_.B = B;
B = min([round(.5*NumberOfDraws),500]); options_.B = B;
end
try delete([MhDirectoryName '/' M_.fname '_irf_dsge*.mat'])
try delete([MhDirectoryName filesep M_.fname '_irf_dsge*.mat'])
catch disp('No _IRFs (dsge) files to be deleted!')
end
try delete([MhDirectoryName '/' M_.fname '_irf_bvardsge*.mat'])
try delete([MhDirectoryName filesep M_.fname '_irf_bvardsge*.mat'])
catch disp('No _IRFs (bvar-dsge) files to be deleted!')
end
irun = 0;
@ -105,23 +106,16 @@ irun2 = 0;
NumberOfIRFfiles_dsge = 1;
NumberOfIRFfiles_dsgevar = 1;
ifil2 = 1;
if strcmpi(type,'posterior')
h = waitbar(0,'Bayesian (posterior) IRFs...');
elseif strcmpi(type,'gsa')
h = waitbar(0,'GSA (prior) IRFs...');
else
h = waitbar(0,'Bayesian (prior) IRFs...');
end
% Create arrays
if B <= MAX_nruns
stock_param = zeros(B, npar);
stock_param = zeros(B, npar);
else
stock_param = zeros(MAX_nruns, npar);
stock_param = zeros(MAX_nruns, npar);
end
if B >= MAX_nirfs_dsge
stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,MAX_nirfs_dsge);
stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,MAX_nirfs_dsge);
else
stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,B);
stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,B);
end
if MAX_nirfs_dsgevar
if B >= MAX_nirfs_dsgevar
@ -142,149 +136,99 @@ if MAX_nirfs_dsgevar
end
Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs);
end
%%%%%%%%% START the FIRST BLOCK of CODE EXECUTED in PARALLEL! %%%%%%%%%
%
% This portion of code is execute in parallel by PosteriorIRF_core1.m
% function.
b = 0;
nosaddle = 0;
while b<=B
localVars=[];
% Save the local variables.
localVars.IRUN = IRUN;
localVars.irun = irun;
localVars.irun2=irun2;
localVars.nosaddle=nosaddle;
% It is necessary to rename 'type' to avoid conflict with
% a native matlab funtion.
localVars.typee=type;
if strcmpi(type,'posterior'),
while b<=B
b = b + 1;
irun = irun+1;
irun2 = irun2+1;
if ~strcmpi(type,'gsa')
deep = GetOneDraw(type);
else
deep = x(b,:);
end
stock_param(irun2,:) = deep;
set_parameters(deep);
[dr,info] = resol(oo_.steady_state,0);
if info(1)
nosaddle = nosaddle + 1;
b = b - 1;
irun = irun-1;
irun2 = irun2-1;
if info(1) == 1
errordef = 'Static variables are not uniquely defined';
elseif info(1) == 2
errordef = 'Dll problem';
elseif info(1) == 3
errordef = 'No stable trajectory';
elseif info(1) == 4
errordef = 'Indeterminacy';
elseif info(1) == 5
errordef = 'Rank condition is not satisfied';
end
disp(['PosteriorIRF :: Dynare is unable to solve the model (' errordef ')'])
continue
end
SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord) = M_.Sigma_e+1e-14*eye(M_.exo_nbr);
SS = transpose(chol(SS));
for i = 1:M_.exo_nbr
if SS(i,i) > 1e-13
y=irf(dr,SS(M_.exo_names_orig_ord,i), options_.irf, options_.drop,options_.replic,options_.order);
if options_.relative_irf
y = 100*y/cs(i,i);
end
for j = 1:nvar
if max(y(IndxVariables(j),:)) - min(y(IndxVariables(j),:)) > 1e-12
stock_irf_dsge(:,j,i,irun) = transpose(y(IndxVariables(j),:));
end
end
end
end
if MAX_nirfs_dsgevar
IRUN = IRUN+1;
%tmp_dsgevar = zeros(options_.irf,nvobs*M_.exo_nbr);
[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)));
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);
% draw from the conditional posterior of PHI
PHI_draw = rand_matrix_normal(NumberOfParametersPerEquation,nvobs, PHI, ...
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
% $$$ if options_.noconstant
mu = zeros(1,nvobs);
% $$$ else
% $$$ AA = eye(nvobs);
% $$$ for lag=1:NumberOfLags
% $$$ AA = AA-PHI_draw((lag-1)*nvobs+1:lag*nvobs,:);
% $$$ end
% $$$ mu = transpose(AA\transpose(PHI_draw(end,:)));
% $$$ end
% Get rotation
if dsge_prior_weight > 0
Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e);
A0 = Atheta(bayestopt_.mfys,:);
[OMEGAstar,SIGMAtr] = qr2(A0');
end
SIGMAu_chol = chol(SIGMAu_draw)';
SIGMAtrOMEGA = SIGMAu_chol*OMEGAstar';
PHIpower = eye(NumberOfLagsTimesNvobs);
irfs = zeros (options_.irf,nvobs*M_.exo_nbr);
tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
irfs(1,:) = tmp3(:)';
for t = 2:options_.irf
PHIpower = Companion_matrix*PHIpower;
tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
irfs(t,:) = tmp3(:)'+kron(ones(1,M_.exo_nbr),mu);
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
tmp_dsgevar(:,j) = (irfs(:,j));
end
end
if IRUN < MAX_nirfs_dsgevar
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,nvobs,M_.exo_nbr);
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;'];,
eval(['save ' instr]);
NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
IRUN =0;
stock_irf_dsgevar = zeros(options_.irf,nvobs,M_.exo_nbr,MAX_nirfs_dsgevar);
end
end
if irun == MAX_nirfs_dsge | irun == B | b == B
if b == B
stock_irf_dsge = stock_irf_dsge(:,:,:,1:irun);
if MAX_nirfs_dsgevar & (b == B | IRUN == B)
stock_irf_bvardsge = stock_irf_bvardsge(:,:,:,1:IRUN);
instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' ...
int2str(NumberOfIRFfiles_dsgevar) ' stock_irf_bvardsge;'];,
eval(['save ' instr]);
NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
irun = 0;
end
end
save([MhDirectoryName '/' M_.fname '_irf_dsge' int2str(NumberOfIRFfiles_dsge) '.mat'],'stock_irf_dsge');
NumberOfIRFfiles_dsge = NumberOfIRFfiles_dsge+1;
irun = 0;
end
if irun2 == MAX_nruns | b == B
if b == B
stock_param = stock_param(1:irun2,:);
end
stock = stock_param;
save([MhDirectoryName '/' M_.fname '_param_irf' int2str(ifil2) '.mat'],'stock');
ifil2 = ifil2 + 1;
irun2 = 0;
end
waitbar(b/B,h);
x(b,:) = GetOneDraw(type);
end
end
if ~strcmpi(type,'prior'),
localVars.x=x;
end
b=0;
localVars.nvar=nvar;
localVars.IndxVariables=IndxVariables;
localVars.MAX_nirfs_dsgevar=MAX_nirfs_dsgevar;
localVars.MAX_nirfs_dsge=MAX_nirfs_dsge;
localVars.MAX_nruns=MAX_nruns;
localVars.NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge;
localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar;
localVars.ifil2=ifil2;
if isnumeric(options_.parallel),% | isunix, % for the moment exclude unix platform from parallel implementation
[fout] = PosteriorIRF_core1(localVars,1,B,0);
else
[nCPU, totCPU, nBlockPerCPU] = distributeJobs(options_.parallel, 1, B);
for j=1:totCPU-1,
nfiles = ceil(nBlockPerCPU(j)/MAX_nirfs_dsge);
NumberOfIRFfiles_dsge(j+1) =NumberOfIRFfiles_dsge(j)+nfiles;
nfiles = ceil(nBlockPerCPU(j)/MAX_nirfs_dsgevar);
NumberOfIRFfiles_dsgevar(j+1) =NumberOfIRFfiles_dsgevar(j)+nfiles;
nfiles = ceil(nBlockPerCPU(j)/MAX_nruns);
ifil2(j+1) =ifil2(j)+nfiles;
end
localVars.NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge;
localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar;
localVars.ifil2=ifil2;
globalVars = struct('M_',M_, ...
'options_', options_, ...
'bayestopt_', bayestopt_, ...
'estim_params_', estim_params_, ...
'oo_', oo_);
% which files have to be copied to run remotely
NamFileInput(1,:) = {'',[M_.fname '_static.m']};
NamFileInput(2,:) = {'',[M_.fname '_dynamic.m']};
if options_.steadystate_flag,
NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '_steadystate.m']};
end
[fout] = masterParallel(options_.parallel, 1, B,NamFileInput,'PosteriorIRF_core1', localVars, globalVars, options_.parallel_info);
end
% END parallel code!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nosaddle
disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))])
end
close(h);
disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))])
end
% if isnumeric(options_.parallel)
% close(h);
% end
ReshapeMatFiles('irf_dsge')
if MAX_nirfs_dsgevar
@ -292,13 +236,13 @@ if MAX_nirfs_dsgevar
end
if strcmpi(type,'gsa')
return
return
end
IRF_DSGEs = dir([MhDirectoryName '/' M_.fname '_IRF_DSGEs*.mat']);
IRF_DSGEs = dir([MhDirectoryName filesep M_.fname '_IRF_DSGEs*.mat']);
NumberOfIRFfiles_dsge = length(IRF_DSGEs);
IRF_BVARDSGEs = dir([MhDirectoryName '/' M_.fname '_IRF_BVARDSGEs*.mat']);
IRF_BVARDSGEs = dir([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs*.mat']);
NumberOfIRFfiles_dsgevar = length(IRF_BVARDSGEs);
@ -310,10 +254,10 @@ DistribIRF = zeros(options_.irf,9,nvar,M_.exo_nbr);
HPDIRF = zeros(options_.irf,2,nvar,M_.exo_nbr);
if options_.TeX
varlist_TeX = [];
for i=1:nvar
varlist_TeX = strvcat(varlist_TeX,M_.endo_names_tex(IndxVariables(i),:));
end
varlist_TeX = [];
for i=1:nvar
varlist_TeX = strvcat(varlist_TeX,M_.endo_names_tex(IndxVariables(i),:));
end
end
fprintf('MH: Posterior (dsge) IRFs...\n');
@ -321,31 +265,32 @@ tit(M_.exo_names_orig_ord,:) = M_.exo_names;
kdx = 0;
for file = 1:NumberOfIRFfiles_dsge
load([MhDirectoryName '/' M_.fname '_IRF_DSGEs' int2str(file) '.mat']);
for i = 1:M_.exo_nbr
for j = 1:nvar
for k = 1:size(STOCK_IRF_DSGE,1)
kk = k+kdx;
[MeanIRF(kk,j,i),MedianIRF(kk,j,i),VarIRF(kk,j,i),HPDIRF(kk,:,j,i),...
DistribIRF(kk,:,j,i)] = posterior_moments(squeeze(STOCK_IRF_DSGE(k,j,i,:)),0,options_.mh_conf_sig);
end
load([MhDirectoryName filesep M_.fname '_IRF_DSGEs' int2str(file) '.mat']);
for i = 1:M_.exo_nbr
for j = 1:nvar
for k = 1:size(STOCK_IRF_DSGE,1)
kk = k+kdx;
[MeanIRF(kk,j,i),MedianIRF(kk,j,i),VarIRF(kk,j,i),HPDIRF(kk,:,j,i),...
DistribIRF(kk,:,j,i)] = posterior_moments(squeeze(STOCK_IRF_DSGE(k,j,i,:)),0,options_.mh_conf_sig);
end
end
kdx = kdx + size(STOCK_IRF_DSGE,1);
end
kdx = kdx + size(STOCK_IRF_DSGE,1);
end
clear STOCK_IRF_DSGE;
for i = 1:M_.exo_nbr
for j = 1:nvar
name = [deblank(M_.endo_names(IndxVariables(j),:)) '_' deblank(tit(i,:))];
eval(['oo_.PosteriorIRF.dsge.Mean.' name ' = MeanIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.dsge.Median.' name ' = MedianIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.dsge.Var.' name ' = VarIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.dsge.Distribution.' name ' = DistribIRF(:,:,j,i);']);
eval(['oo_.PosteriorIRF.dsge.HPDinf.' name ' = HPDIRF(:,1,j,i);']);
eval(['oo_.PosteriorIRF.dsge.HPDsup.' name ' = HPDIRF(:,2,j,i);']);
end
for j = 1:nvar
name = [deblank(M_.endo_names(IndxVariables(j),:)) '_' deblank(tit(i,:))];
eval(['oo_.PosteriorIRF.dsge.Mean.' name ' = MeanIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.dsge.Median.' name ' = MedianIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.dsge.Var.' name ' = VarIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.dsge.Distribution.' name ' = DistribIRF(:,:,j,i);']);
eval(['oo_.PosteriorIRF.dsge.HPDinf.' name ' = HPDIRF(:,1,j,i);']);
eval(['oo_.PosteriorIRF.dsge.HPDsup.' name ' = HPDIRF(:,2,j,i);']);
end
end
@ -359,7 +304,7 @@ if MAX_nirfs_dsgevar
tit(M_.exo_names_orig_ord,:) = M_.exo_names;
kdx = 0;
for file = 1:NumberOfIRFfiles_dsgevar
load([MhDirectoryName '/' M_.fname '_IRF_BVARDSGEs' int2str(file) '.mat']);
load([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs' int2str(file) '.mat']);
for i = 1:M_.exo_nbr
for j = 1:nvar
for k = 1:size(STOCK_IRF_BVARDSGE,1)
@ -386,114 +331,101 @@ if MAX_nirfs_dsgevar
end
end
%%
%% Finally I build the plots.
%% Finally I build the plots.
%%
%%%%%%%%% START the SECOND BLOCK of CODE EXECUTED in PARALLEL! %%%%%%%%%
%
% This portion of code is execute in parallel by PosteriorIRF_core2.m
% function.
% Save the local variables.
localVars=[];
Check=options_.TeX
if (Check)
localVars.varlist_TeX=varlist_TeX;
end
localVars.nvar=nvar;
localVars.MeanIRF=MeanIRF;
localVars.tit=tit;
localVars.nn=nn;
localVars.MAX_nirfs_dsgevar=MAX_nirfs_dsgevar;
localVars.HPDIRF=HPDIRF;
localVars.varlist=varlist;
localVars.MaxNumberOfPlotPerFigure=MaxNumberOfPlotPerFigure;
%%% The files .TeX are genereted in sequential way!
if options_.TeX
fidTeX = fopen([DirectoryName '/' M_.fname '_BayesianIRF.TeX'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by PosteriorIRF.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
titTeX(M_.exo_names_orig_ord,:) = M_.exo_names_tex;
end
%%
subplotnum = 0;
for i=1:M_.exo_nbr
NAMES = [];
if options_.TeX
fidTeX = fopen([DirectoryName filesep M_.fname '_BayesianIRF.TeX'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by PosteriorIRF.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
titTeX(M_.exo_names_orig_ord,:) = M_.exo_names_tex;
for i=1:M_.exo_nbr
NAMES = [];
TEXNAMES = [];
end
figunumber = 0;
for j=1:nvar
if max(abs(MeanIRF(:,j,i))) > 10^(-6)
subplotnum = subplotnum+1;
if options_.nograph
if subplotnum == 1 & options_.relative_irf
hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)],'Visible','off');
elseif subplotnum == 1 & ~options_.relative_irf
hh = figure('Name',['Orthogonalized shock to ' tit(i,:)],'Visible','off');
end
else
if subplotnum == 1 & options_.relative_irf
hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)]);
elseif subplotnum == 1 & ~options_.relative_irf
hh = figure('Name',['Orthogonalized shock to ' tit(i,:)]);
end
end
set(0,'CurrentFigure',hh)
subplot(nn,nn,subplotnum);
if ~MAX_nirfs_dsgevar
h1 = area(1:options_.irf,HPDIRF(:,2,j,i));
set(h1,'FaceColor',[.9 .9 .9]);
set(h1,'BaseValue',min(HPDIRF(:,1,j,i)));
hold on
h2 = area(1:options_.irf,HPDIRF(:,1,j,i),'FaceColor',[1 1 1],'BaseValue',min(HPDIRF(:,1,j,i)));
set(h2,'FaceColor',[1 1 1]);
set(h2,'BaseValue',min(HPDIRF(:,1,j,i)));
plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',3)
% plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
box on
axis tight
xlim([1 options_.irf]);
hold off
else
h1 = area(1:options_.irf,HPDIRF(:,2,j,i));
set(h1,'FaceColor',[.9 .9 .9]);
set(h1,'BaseValue',min([min(HPDIRF(:,1,j,i)),min(HPDIRFdsgevar(:,1,j,i))]));
hold on;
h2 = area(1:options_.irf,HPDIRF(:,1,j,i));
set(h2,'FaceColor',[1 1 1]);
set(h2,'BaseValue',min([min(HPDIRF(:,1,j,i)),min(HPDIRFdsgevar(:,1,j,i))]));
plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',3)
% plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
plot(1:options_.irf,MeanIRFdsgevar(:,j,i),'--k','linewidth',2)
plot(1:options_.irf,HPDIRFdsgevar(:,1,j,i),'--k','linewidth',1)
plot(1:options_.irf,HPDIRFdsgevar(:,2,j,i),'--k','linewidth',1)
box on
axis tight
xlim([1 options_.irf]);
hold off
end
for j=1:nvar
if max(abs(MeanIRF(:,j,i))) > 10^(-6)
name = deblank(varlist(j,:));
NAMES = strvcat(NAMES,name);
if options_.TeX
texname = deblank(varlist_TeX(j,:));
TEXNAMES = strvcat(TEXNAMES,['$' texname '$']);
texname = deblank(varlist_TeX(j,:));
TEXNAMES = strvcat(TEXNAMES,['$' texname '$']);
end
end
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:size(TEXNAMES,1)
fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{%s}\n'],deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:)));
end
title(name,'Interpreter','none')
end
if subplotnum == MaxNumberOfPlotPerFigure | (j == nvar & subplotnum> 0)
figunumber = figunumber+1;
set(hh,'visible','on')
eval(['print -depsc2 ' DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.eps']);
if ~exist('OCTAVE_VERSION')
eval(['print -dpdf ' DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)]);
saveas(hh,[DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.fig']);
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_Bayesian_IRF_%s}\n',M_.fname,deblank(tit(i,:)));
if options_.relative_irf
fprintf(fidTeX,['\\caption{Bayesian relative IRF.}']);
else
fprintf(fidTeX,'\\caption{Bayesian IRF.}');
end
set(hh,'visible','off')
if options_.nograph, close(hh), end
if options_.TeX
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:size(TEXNAMES,1)
fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{%s}\n'],deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:)));
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_Bayesian_IRF_%s}\n',M_.fname,deblank(tit(i,:)));
if options_.relative_irf
fprintf(fidTeX,['\\caption{Bayesian relative IRF.}']);
else
fprintf(fidTeX,'\\caption{Bayesian IRF.}');
end
fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s}\n',deblank(tit(i,:)));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
subplotnum = 0;
end
end% loop over selected endo_var
end% loop over exo_var
%%
if options_.TeX
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s}\n',deblank(tit(i,:)));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
end
fprintf('MH: Posterior IRFs, done!\n');
% The others format in parallel by PosteriorIRF_core2!
if isnumeric(options_.parallel) || (M_.exo_nbr*ceil(size(varlist,1)/MaxNumberOfPlotPerFigure))<8,% | isunix, % for the moment exclude unix platform from parallel implementation
[fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0);
else
globalVars = struct('M_',M_, ...
'options_', options_);
[fout] = masterParallel(options_.parallel, 1, M_.exo_nbr,NamFileInput,'PosteriorIRF_core2', localVars, globalVars, options_.parallel_info);
end
% END parallel code!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('MH: Posterior IRFs, done!\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

278
matlab/PosteriorIRF_core1.m Normal file
View File

@ -0,0 +1,278 @@
function myoutput=PosteriorIRF_core1(myinputs,fpar,npar,whoiam, ThisMatlab)
% Perfome in parallel a portion of PosteriorIRF
%
% INPUTS
%
% ...
%
% OUTPUTS
% ...
%
% SPECIAL REQUIREMENTS
% None
% Copyright (C) 2006-2008 Dynare Team
%
% This file is part of Dynare.
%
% 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/>.
global options_ estim_params_ oo_ M_ bayestopt_
if nargin<4,
whoiam=0;
end
struct2local(myinputs);
MhDirectoryName = CheckPath('metropolis');
RemoteFlag = 0;
if whoiam
waitbarString = ['Please wait... Bayesian (posterior) IRFs computing. (' int2str(fpar) 'of' int2str(npar) ')...'];
if Parallel(ThisMatlab).Local,
waitbarTitle=['Local '];
else
waitbarTitle=[Parallel(ThisMatlab).PcName];
RemoteFlag =1;
end
fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo);
else
if exist('OCTAVE_VERSION')
diary off;
printf('\n')
else
if strcmpi(typee,'posterior')
h = waitbar(0,'Bayesian (posterior) IRFs...');
set(h,'Name','Bayesian (posterior) IRFs.');
elseif strcmpi(typee,'gsa')
h = waitbar(0,'GSA (prior) IRFs...');
else
h = waitbar(0,'Bayesian (prior) IRFs...');
end
end
end
OutputFileName_bvardsge = {};
OutputFileName_dsge = {};
OutputFileName_param = {};
fpar0=fpar;
fpar = fpar-1;
if whoiam
ifil2=ifil2(whoiam);
NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge(whoiam);
NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar(whoiam);
end
while fpar<=npar
fpar = fpar + 1;
irun = irun+1;
irun2 = irun2+1;
if strcmpi(typee,'prior')
deep = GetOneDraw(typee);
else
deep = x(fpar,:);
end
stock_param(irun2,:) = deep;
set_parameters(deep);
[dr,info] = resol(oo_.steady_state,0);
if info(1)
nosaddle = nosaddle + 1;
fpar = fpar - 1;
irun = irun-1;
irun2 = irun2-1;
if info(1) == 1
errordef = 'Static variables are not uniquely defined';
elseif info(1) == 2
errordef = 'Dll problem';
elseif info(1) == 3
errordef = 'No stable trajectory';
elseif info(1) == 4
errordef = 'Indeterminacy';
elseif info(1) == 5
errordef = 'Rank condition is not satisfied';
end
disp(['PosteriorIRF :: Dynare is unable to solve the model (' errordef ')'])
continue
end
SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord) = M_.Sigma_e+1e-14*eye(M_.exo_nbr);
SS = transpose(chol(SS));
for i = 1:M_.exo_nbr
if SS(i,i) > 1e-13
y=irf(dr,SS(M_.exo_names_orig_ord,i), options_.irf, options_.drop,options_.replic,options_.order);
if options_.relative_irf
y = 100*y/cs(i,i);
end
for j = 1:nvar
if max(y(IndxVariables(j),:)) - min(y(IndxVariables(j),:)) > 1e-12
stock_irf_dsge(:,j,i,irun) = transpose(y(IndxVariables(j),:));
end
end
end
end
if MAX_nirfs_dsgevar
IRUN = IRUN+1;
%tmp_dsgevar = zeros(options_.irf,nvobs*M_.exo_nbr);
[fval,cost_flag,info,PHI,SIGMAu,iXX] = DsgeVarLikelihood(deep(fpar,:)',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)));
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);
% draw from the conditional posterior of PHI
PHI_draw = rand_matrix_normal(NumberOfParametersPerEquation,nvobs, PHI, ...
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
% $$$ if options_.noconstant
mu = zeros(1,nvobs);
% $$$ else
% $$$ AA = eye(nvobs);
% $$$ for lag=1:NumberOfLags
% $$$ AA = AA-PHI_draw((lag-1)*nvobs+1:lag*nvobs,:);
% $$$ end
% $$$ mu = transpose(AA\transpose(PHI_draw(end,:)));
% $$$ end
% Get rotation
if dsge_prior_weight > 0
Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e);
A0 = Atheta(bayestopt_.mfys,:);
[OMEGAstar,SIGMAtr] = qr2(A0');
end
SIGMAu_chol = chol(SIGMAu_draw)';
SIGMAtrOMEGA = SIGMAu_chol*OMEGAstar';
PHIpower = eye(NumberOfLagsTimesNvobs);
irfs = zeros (options_.irf,nvobs*M_.exo_nbr);
tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
irfs(1,:) = tmp3(:)';
for t = 2:options_.irf
PHIpower = Companion_matrix*PHIpower;
tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
irfs(t,:) = tmp3(:)'+kron(ones(1,M_.exo_nbr),mu);
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
tmp_dsgevar(:,j) = (irfs(:,j));
end
end
if IRUN < MAX_nirfs_dsgevar
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,nvobs,M_.exo_nbr);
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;'];,
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;
IRUN =0;
stock_irf_dsgevar = zeros(options_.irf,nvobs,M_.exo_nbr,MAX_nirfs_dsgevar);
end
end
if irun == MAX_nirfs_dsge | irun == npar | fpar == npar
if fpar == npar
stock_irf_dsge = stock_irf_dsge(:,:,:,1:irun);
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) ' stock_irf_bvardsge;'];,
eval(['save ' instr]);
NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
if RemoteFlag==1,
OutputFileName_bvardsge = [OutputFileName_bvardsge; {[MhDirectoryName filesep], [M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) '.mat']}];
end
irun = 0;
end
end
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
NumberOfIRFfiles_dsge = NumberOfIRFfiles_dsge+1;
irun = 0;
end
if irun2 == MAX_nruns | fpar == npar
if fpar == npar
stock_param = stock_param(1:irun2,:);
end
stock = stock_param;
save([MhDirectoryName '/' M_.fname '_param_irf' int2str(ifil2) '.mat'],'stock');
if RemoteFlag==1,
OutputFileName_param = [OutputFileName_param; {[MhDirectoryName filesep], [M_.fname '_param_irf' int2str(ifil2) '.mat']}];
end
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);
end
if mod(fpar,10)==0 & whoiam,
fprintf('Done! \n');
waitbarString = [ 'Subdraw ' int2str(fpar) '/' int2str(npar) ' done.'];
fMessageStatus((fpar-fpar0)/(npar-fpar0),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo)
end
end
if whoiam==0
if nosaddle
disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(npar+nosaddle))])
end
if exist('h')
close(h);
end
if exist('OCTAVE_VERSION')
printf('\n');
diary on;
end
end
% Copy the rusults of computation on the call machine (specifically in the
% dyrectory on call machine that contain the model).
myoutput.OutputFileName = [OutputFileName_dsge;
OutputFileName_param;
OutputFileName_bvardsge];
% MhDirectoryName=TempPath;

121
matlab/PosteriorIRF_core2.m Normal file
View File

@ -0,0 +1,121 @@
function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab)
global options_ M_
if nargin<4,
whoiam=0;
end
struct2local(myinputs);
% To save the figures where the function is computed!
DirectoryName = CheckPath('Output');
RemoteFlag = 0;
if whoiam,
waitbarString = ['Please wait... PosteriorIRF Plots (exog. shocks ' int2str(fpar) 'of' int2str(npar) ')...'];
if Parallel(ThisMatlab).Local,
waitbarTitle=['Local '];
else
waitbarTitle=[Parallel(ThisMatlab).PcName];
RemoteFlag = 1;
end
fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo);
end
OutputFileName={};
subplotnum = 0;
for i=fpar:npar,
NAMES = [];
figunumber = 0;
for j=1:nvar
if max(abs(MeanIRF(:,j,i))) > 10^(-6)
subplotnum = subplotnum+1;
if options_.nograph
if subplotnum == 1 & options_.relative_irf
hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)],'Visible','off');
elseif subplotnum == 1 & ~options_.relative_irf
hh = figure('Name',['Orthogonalized shock to ' tit(i,:)],'Visible','off');
end
else
if subplotnum == 1 & options_.relative_irf
hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)]);
elseif subplotnum == 1 & ~options_.relative_irf
hh = figure('Name',['Orthogonalized shock to ' tit(i,:)]);
end
end
set(0,'CurrentFigure',hh)
subplot(nn,nn,subplotnum);
if ~MAX_nirfs_dsgevar
h1 = area(1:options_.irf,HPDIRF(:,2,j,i));
set(h1,'FaceColor',[.9 .9 .9]);
set(h1,'BaseValue',min(HPDIRF(:,1,j,i)));
hold on
h2 = area(1:options_.irf,HPDIRF(:,1,j,i),'FaceColor',[1 1 1],'BaseValue',min(HPDIRF(:,1,j,i)));
set(h2,'FaceColor',[1 1 1]);
set(h2,'BaseValue',min(HPDIRF(:,1,j,i)));
plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',3)
% plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
box on
axis tight
xlim([1 options_.irf]);
hold off
else
h1 = area(1:options_.irf,HPDIRF(:,2,j,i));
set(h1,'FaceColor',[.9 .9 .9]);
set(h1,'BaseValue',min([min(HPDIRF(:,1,j,i)),min(HPDIRFdsgevar(:,1,j,i))]));
hold on;
h2 = area(1:options_.irf,HPDIRF(:,1,j,i));
set(h2,'FaceColor',[1 1 1]);
set(h2,'BaseValue',min([min(HPDIRF(:,1,j,i)),min(HPDIRFdsgevar(:,1,j,i))]));
plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',3)
% plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
plot(1:options_.irf,MeanIRFdsgevar(:,j,i),'--k','linewidth',2)
plot(1:options_.irf,HPDIRFdsgevar(:,1,j,i),'--k','linewidth',1)
plot(1:options_.irf,HPDIRFdsgevar(:,2,j,i),'--k','linewidth',1)
box on
axis tight
xlim([1 options_.irf]);
hold off
end
name = deblank(varlist(j,:));
NAMES = strvcat(NAMES,name);
title(name,'Interpreter','none')
end
% TempPath=DirectoryName;
% DirectoryNamePar='C:\dynare_calcs\ModelTest\ls2003\metropolis';
% DirectoryName=DirectoryNamePar;
if subplotnum == MaxNumberOfPlotPerFigure | (j == nvar & subplotnum> 0)
figunumber = figunumber+1;
set(hh,'visible','on')
eval(['print -depsc2 ' DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.eps']);
if ~exist('OCTAVE_VERSION')
eval(['print -dpdf ' DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)]);
saveas(hh,[DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.fig']);
end
if RemoteFlag==1,
OutputFileName = [OutputFileName; {[DirectoryName,filesep], [M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.*']}];
end
set(hh,'visible','off')
if options_.nograph, close(hh), end
subplotnum = 0;
end
end% loop over selected endo_var
if whoiam,
fprintf('Done! \n');
waitbarString = [ 'Exog. shocks ' int2str(i) '/' int2str(npar) ' done.'];
fMessageStatus((i-fpar+1)/(npar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo)
end
end% loop over exo_var
% DirectoryName=TempPath;
myoutput.OutputFileName = OutputFileName;

12
matlab/ReshapeMatFiles.m
View File

@ -45,20 +45,20 @@ function ReshapeMatFiles(type, type2)
global M_ options_
if nargin==1,
MhDirectoryName = [ CheckPath('metropolis') '/' ];
MhDirectoryName = [ CheckPath('metropolis') filesep ];
else
if strcmpi(type2,'posterior')
MhDirectoryName = [CheckPath('metropolis') '/' ];
MhDirectoryName = [CheckPath('metropolis') filesep ];
elseif strcmpi(type2,'gsa')
if options_.opt_gsa.morris==1,
MhDirectoryName = [CheckPath('GSA\SCREEN') '/' ];
MhDirectoryName = [CheckPath('GSA\SCREEN') filesep ];
elseif options_.opt_gsa.morris==2,
MhDirectoryName = [CheckPath('GSA\IDENTIF') '/' ];
MhDirectoryName = [CheckPath('GSA\IDENTIF') filesep ];
else
MhDirectoryName = [CheckPath('GSA') '/' ];
MhDirectoryName = [CheckPath('GSA') filesep ];
end
else
MhDirectoryName = [CheckPath('prior') '/' ];
MhDirectoryName = [CheckPath('prior') filesep ];
end
end
switch type

40
matlab/Tracing.m Normal file
View File

@ -0,0 +1,40 @@
function [] = Tracing()
% DESCRIPTION
% This function is used to test the correct execution of a matlab section
% on remote machine.
%
% If no error happen the function simply create a file.
% INPUTS
% ...
%
% OUTPUTS
% ...
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2005-2009 Dynare Team
%
% This file is part of Dynare.
%
% 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/>.
fid = fopen('Tracing.txt','w+');
fclose (fid);
exit

14
matlab/closeSlave.m Normal file
View File

@ -0,0 +1,14 @@
function closeSlave(Parallel),
% In a parallelc context, this utility closes all remote matlab instances
% called by masteParallelMan (which leaves open remote matlab instances)
delete( 'slaveParallel_input*.mat');
for indPC=1:length(Parallel),
if Parallel(indPC).Local==0,
if isunix || (~matlab_ver_less_than('7.4') && ismac),
system(['ssh ',Parallel(indPC).user,'@',Parallel(indPC).PcName,' rm -fr ',Parallel(indPC).RemoteFolder,'/slaveParallel_input*.mat']);
else
mydelete('slaveParallel_input*.mat',['\\',Parallel(indPC).PcName,'\',Parallel(indPC).RemoteDrive,'$\',Parallel(indPC).RemoteFolder,'\']);
end
end
end

37
matlab/distributeJobs.m Normal file
View File

@ -0,0 +1,37 @@
function [nCPU, totCPU, nBlockPerCPU] = distributeJobs(Parallel, fBlock, nBlock)
% Determine the total number of available CPUs, and the number of threads to run on each CPU
% Copyright (C) 2009 Dynare Team
%
% This file is part of Dynare.
%
% 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/>.
totCPU=0;
for j=1:length(Parallel),
nCPU(j)=length(Parallel(j).NumCPU);
totCPU=totCPU+nCPU(j);
end
nCPU=cumsum(nCPU);
offset0 = fBlock-1;
if (nBlock-offset0)>totCPU,
diff = mod((nBlock-offset0),totCPU);
nBlockPerCPU(1:diff) = ceil((nBlock-offset0)/totCPU);
nBlockPerCPU(diff+1:totCPU) = floor((nBlock-offset0)/totCPU);
else
nBlockPerCPU(1:nBlock-offset0)=1;
totCPU = nBlock-offset0;
end

1
matlab/global_initialization.m
View File

@ -191,6 +191,7 @@ options_.mh_autocorrelation_function_size = 30;
options_.plot_priors = 1;
options_.cova_compute = 1;
options_.parallel = 0;
options_.parallel_info.leaveSlaveOpen = 0;
options_.number_of_grid_points_for_kde = 2^9;
quarter = 1;
years = [1 2 3 4 5 10 20 30 40 50];

2
matlab/independent_metropolis_hastings.m
View File

@ -95,7 +95,7 @@ else
% from where to get back results
% NamFileOutput(1,:) = {[M_.dname,'/metropolis/'],'*.*'};
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'independent_metropolis_hastings_core', localVars, globalVars);
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'independent_metropolis_hastings_core', localVars, globalVars, options_.parallel_info);
for j=1:totCPU,
offset = sum(nBlockPerCPU(1:j-1))+fblck-1;
record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)));

48
matlab/masterParallel.m
View File

@ -1,4 +1,4 @@
function [fOutVar,nBlockPerCPU, totCPU] = masterParallel(Parallel,fBlock,nBlock,NamFileInput,fname,fInputVar,fGlobalVar)
function [fOutVar,nBlockPerCPU, totCPU] = masterParallel(Parallel,fBlock,nBlock,NamFileInput,fname,fInputVar,fGlobalVar,Parallel_info)
% Top-level function called on the master computer when parallelizing a task.
%
% The number of parallelized threads will be equal to (nBlock-fBlock+1).
@ -46,9 +46,13 @@ function [fOutVar,nBlockPerCPU, totCPU] = masterParallel(Parallel,fBlock,nBlock,
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
totCPU=0;
% Determine my hostname and my working directory
if ~isempty(Parallel_info)
if Parallel_info.leaveSlaveOpen,
[fOutVar,nBlockPerCPU, totCPU] = masterParallelMan(Parallel,fBlock,nBlock,NamFileInput,fname,fInputVar,fGlobalVar);
return
end
end
DyMo=pwd;
fInputVar.DyMo=DyMo;
if isunix || (~matlab_ver_less_than('7.4') && ismac) ,
@ -69,21 +73,24 @@ end
save([fname,'_input.mat'],'Parallel','-append')
% Determine the total number of available CPUs, and the number of threads to run on each CPU
for j=1:length(Parallel),
nCPU(j)=length(Parallel(j).NumCPU);
totCPU=totCPU+nCPU(j);
end
nCPU=cumsum(nCPU);
[nCPU, totCPU, nBlockPerCPU] = distributeJobs(Parallel, fBlock, nBlock);
offset0 = fBlock-1;
if (nBlock-offset0)>totCPU,
diff = mod((nBlock-offset0),totCPU);
nBlockPerCPU(1:diff) = ceil((nBlock-offset0)/totCPU);
nBlockPerCPU(diff+1:totCPU) = floor((nBlock-offset0)/totCPU);
else
nBlockPerCPU(1:nBlock-offset0)=1;
totCPU = nBlock-offset0;
end
% for j=1:length(Parallel),
% nCPU(j)=length(Parallel(j).NumCPU);
% totCPU=totCPU+nCPU(j);
% end
%
% nCPU=cumsum(nCPU);
% offset0 = fBlock-1;
% if (nBlock-offset0)>totCPU,
% diff = mod((nBlock-offset0),totCPU);
% nBlockPerCPU(1:diff) = ceil((nBlock-offset0)/totCPU);
% nBlockPerCPU(diff+1:totCPU) = floor((nBlock-offset0)/totCPU);
% else
% nBlockPerCPU(1:nBlock-offset0)=1;
% totCPU = nBlock-offset0;
% end
% Clean up remnants of previous runs
mydelete(['comp_status_',fname,'*.mat'])
@ -214,6 +221,7 @@ t00=cputime;
hh=NaN(1,nBlock);
if exist('OCTAVE_VERSION'),
diary off;
printf('\n');
else
hfigstatus = figure('name',['Parallel ',fname],...
'MenuBar', 'none', ...
@ -241,7 +249,7 @@ while (1)
load(stax(j).name)
pcerdone(j) = prtfrc;
if exist('OCTAVE_VERSION'),
statusString = [statusString, waitbarString, ', %3.f%% done! '];
statusString = [statusString, int2str(j), ' %3.f%% done! '];
else
status_String{j} = waitbarString;
status_Title{j} = waitbarTitle;
@ -257,6 +265,7 @@ while (1)
else
figure(hfigstatus),
for j=1:length(stax),
try
axes(hstatus(idCPU(j))),
hpat = findobj(hstatus(idCPU(j)),'Type','patch');
if ~isempty(hpat),
@ -265,6 +274,9 @@ while (1)
patch([0 0 pcerdone(j) pcerdone(j)],[0 1 1 0],'r','EdgeColor','r')
end
title([status_Title{j},' - ',status_String{j}]);
catch
end
end
end
if isempty(dir(['P_',fname,'_*End.txt']))

7
matlab/masterParallelMan.m
View File

@ -52,6 +52,7 @@ persistent initialize
if isempty(initialize),
mydelete(['P_slave_*End.txt']);
mydelete(['slaveParallel_input*.mat']);
initialize = 0;
pause(1),
end
@ -101,6 +102,7 @@ mydelete(['P_',fname,'*End.txt']);
% Create a shell script containing the commands to launch the required tasks on the slaves
fid = fopen('ConcurrentCommand1.bat','w+');
for j=1:totCPU,
command1 = ' ';
indPC=min(find(nCPU>=j));
@ -156,7 +158,7 @@ for j=1:totCPU,
command1=['start /B psexec -W ',DyMo, ' -a ',int2str(Parallel(indPC).NumCPU(j-nCPU0)),' -low matlab -nosplash -nodesktop -minimize -r slaveParallel(',int2str(j),',',int2str(indPC),')'];
end
end
else
elseif Parallel(indPC).Local==0,
if isunix || (~matlab_ver_less_than('7.4') && ismac),
% [tempo, RemoteName]=system(['ssh ',Parallel(indPC).user,'@',Parallel(indPC).PcName,' "ifconfig | grep \''inet addr:\''| grep -v \''127.0.0.1\'' | cut -d: -f2 | awk \''{ print $1}\''"']);
[tempo, RemoteName]=system(['ssh ',Parallel(indPC).user,'@',Parallel(indPC).PcName,' "hostname --fqdn"']);
@ -263,6 +265,7 @@ t00=cputime;
hh=NaN(1,nBlock);
if exist('OCTAVE_VERSION'),
diary off;
printf('\n');
else
hfigstatus = figure('name',['Parallel ',fname],...
'MenuBar', 'none', ...
@ -290,7 +293,7 @@ while (1)
load(stax(j).name)
pcerdone(j) = prtfrc;
if exist('OCTAVE_VERSION'),
statusString = [statusString, waitbarString, ', %3.f%% done! '];
statusString = [statusString, int2str(j), ' %3.f%% done! '];
else
status_String{j} = waitbarString;
status_Title{j} = waitbarTitle;

105
matlab/pm3.m
View File

@ -77,77 +77,80 @@ for i = 1:nvar
eval(['oo_.' name3 '.HPDsup.' name ' = HPD(2,:,i);']);
end
%%
%% Finally I build the plots.
%% Finally I build the plots.
%%
% %%%%%%%%% PARALLEL BLOCK % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% %%% The file .TeX! are not saved in parallel.
subplotnum = 0;
if options_.TeX
fidTeX = fopen([M_.dname '/Output/' M_.fname '_' name3 '.TeX'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by Dynare.\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
end
%%
figunumber = 0;
subplotnum = 0;
hh = figure('Name',[tit1 ' ' int2str(figunumber+1)]);
for i=1:nvar
NAMES = [];
if options_.TeX
TEXNAMES = [];
end
if max(abs(Mean(:,i))) > 10^(-6)
subplotnum = subplotnum+1;
set(0,'CurrentFigure',hh)
subplot(nn,nn,subplotnum);
plot([1 n2],[0 0],'-r','linewidth',0.5);
hold on
for k = 1:9
plot(1:n2,squeeze(Distrib(k,:,i)),'-g','linewidth',0.5)
end
plot(1:n2,Mean(:,i),'-k','linewidth',1)
xlim([1 n2]);
hold off
name = deblank(varlist(i,:));
NAMES = strvcat(NAMES,name);
if options_.TeX
for i=1:nvar
NAMES = [];
TEXNAMES = [];
if max(abs(Mean(:,i))) > 10^(-6)
subplotnum = subplotnum+1;
name = deblank(varlist(i,:));
NAMES = strvcat(NAMES,name);
texname = deblank(varlist_TeX(i,:));
TEXNAMES = strvcat(TEXNAMES,['$' texname '$']);
end
title(name,'Interpreter','none')
end
if subplotnum == MaxNumberOfPlotsPerFigure | i == nvar
eval(['print -depsc2 ' M_.dname '/Output/' M_.fname '_' name3 '_' deblank(tit3(i,:)) '.eps' ]);
if ~exist('OCTAVE_VERSION')
eval(['print -dpdf ' M_.dname '/Output/' M_.fname '_' name3 '_' deblank(tit3(i,:))]);
saveas(hh,[M_.dname '/Output/' M_.fname '_' name3 '_' deblank(tit3(i,:)) '.fig']);
end
if options_.nograph, close(hh), end
if options_.TeX
if subplotnum == MaxNumberOfPlotsPerFigure | i == nvar
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:size(TEXNAMES,1)
fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{%s}\n'],deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:)));
end
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,['\\includegraphics[scale=0.5]{%s_' name3 '_%s}\n'],M_.fname,deblank(tit3(i,:)));
if options_.relative_irf
fprintf(fidTeX,['\\caption{' caption '.}']);
else
fprintf(fidTeX,['\\caption{' caption '.}']);
end
fprintf(fidTeX,'\\label{Fig:%s:%s}\n',name3,deblank(tit3(i,:)));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
subplotnum = 0;
figunumber = figunumber+1;
if (i ~= nvar)
hh = figure('Name',[name3 ' ' int2str(figunumber+1)]);
subplotnum = 0;
end
end
end
%%
if options_.TeX
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
end
% Store the variable mandatory for local/remote parallel computing.
localVars=[];
localVars.tit1=tit1;
localVars.nn=nn;
localVars.n2=n2;
localVars.Distrib=Distrib;
localVars.varlist=varlist;
localVars.MaxNumberOfPlotsPerFigure=MaxNumberOfPlotsPerFigure;
localVars.name3=name3;
localVars.tit3=tit3;
localVars.Mean=Mean;
if isnumeric(options_.parallel) || ceil(size(varlist,1)/MaxNumberOfPlotsPerFigure)<4,
fout = pm3_core(localVars,1,nvar,0);
else
globalVars = struct('M_',M_, ...
'options_', options_, ...
'oo_', oo_);
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, 1, nvar, [],'pm3_core', localVars,globalVars, options_.parallel_info);
end
%%%%%%%%% END PARALLEL BLOCK %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(['MH: ' tit1 ', done!\n']);

81
matlab/pm3_core.m Normal file
View File

@ -0,0 +1,81 @@
function myoutput=pm3_core(myinputs,fpar,nvar,whoiam, ThisMatlab)
if nargin<4,
whoiam=0;
end
struct2local(myinputs);
global options_ M_ oo_
if whoiam
waitbarString = ['Parallel plots pm3 ...'];
if Parallel(ThisMatlab).Local,
waitbarTitle=['Local '];
else
waitbarTitle=[Parallel(ThisMatlab).PcName];
end
fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo);
end
figunumber = 0;
subplotnum = 0;
hh = figure('Name',[tit1 ' ' int2str(figunumber+1)]);
RemoteFlag = 0;
if whoiam,
if Parallel(ThisMatlab).Local ==0
RemoteFlag=1;
end
end
OutputFileName = {};
for i=fpar:nvar
NAMES = [];
if max(abs(Mean(:,i))) > 10^(-6)
subplotnum = subplotnum+1;
set(0,'CurrentFigure',hh)
subplot(nn,nn,subplotnum);
plot([1 n2],[0 0],'-r','linewidth',0.5);
hold on
for k = 1:9
plot(1:n2,squeeze(Distrib(k,:,i)),'-g','linewidth',0.5)
end
plot(1:n2,Mean(:,i),'-k','linewidth',1)
xlim([1 n2]);
hold off
name = deblank(varlist(i,:));
NAMES = strvcat(NAMES,name);
title(name,'Interpreter','none')
end
if subplotnum == MaxNumberOfPlotsPerFigure | i == nvar
eval(['print -depsc2 ' M_.dname '/Output/' M_.fname '_' name3 '_' deblank(tit3(i,:)) '.eps' ]);
if ~exist('OCTAVE_VERSION')
eval(['print -dpdf ' M_.dname '/Output/' M_.fname '_' name3 '_' deblank(tit3(i,:))]);
saveas(hh,[M_.dname '/Output/' M_.fname '_' name3 '_' deblank(tit3(i,:)) '.fig']);
end
if RemoteFlag==1,
OutputFileName = [OutputFileName; {[M_.dname, filesep, 'Output',filesep], [M_.fname '_' name3 '_' deblank(tit3(i,:)) '.*']}];
end
if options_.nograph, close(hh), end
subplotnum = 0;
figunumber = figunumber+1;
if (i ~= nvar)
hh = figure('Name',[name3 ' ' int2str(figunumber+1)]);
end
end
if whoiam,
waitbarString = [ 'Variable ' int2str(i) '/' int2str(nvar) ' done.'];
fMessageStatus((i-fpar+1)/(nvar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo)
end
end
myoutput.OutputFileName=OutputFileName;

263
matlab/prior_posterior_statistics.m
View File

@ -1,4 +1,5 @@
function prior_posterior_statistics(type,Y,gend,data_index,missing_value)
% function PosteriorFilterSmootherAndForecast(Y,gend, type)
% Computes posterior filter smoother and forecasts
%
@ -7,13 +8,13 @@ function prior_posterior_statistics(type,Y,gend,data_index,missing_value)
% prior
% gsa
% Y: data
% gend: number of observations
% gend: number of observations
% data_index [cell] 1*smpl cell of column vectors of indices.
% missing_value 1 if missing values, 0 otherwise
%
%
% OUTPUTS
% none
%
%
% SPECIAL REQUIREMENTS
% none
@ -36,6 +37,8 @@ function prior_posterior_statistics(type,Y,gend,data_index,missing_value)
global options_ estim_params_ oo_ M_ bayestopt_
localVars=[];
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
@ -61,8 +64,8 @@ maxlag = M_.maximum_endo_lag;
DirectoryName = CheckPath('metropolis');
load([ DirectoryName '/' M_.fname '_mh_history'])
FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); LastMhFile = TotalNumberOfMhFiles;
FirstLine = record.KeepedDraws.FirstLine;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); LastMhFile = TotalNumberOfMhFiles;
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
clear record;
@ -82,14 +85,14 @@ MAX_ninno = min(B,ceil(MaxNumberOfBytes/(exo_nbr*gend)/8));
MAX_nerro = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)*gend)/8));
if naK
MAX_naK = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)* ...
length(options_.filter_step_ahead)*gend)/8));
length(options_.filter_step_ahead)*gend)/8));
end
if horizon
MAX_nforc1 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/8));
MAX_nforc2 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/ ...
8));
8));
IdObs = bayestopt_.mfys;
end
MAX_momentsno = min(B,ceil(MaxNumberOfBytes/(get_moments_size(options_)*8)));
%%
@ -108,11 +111,6 @@ end
irun = ones(7,1);
ifil = zeros(7,1);
if exist('OCTAVE_VERSION')
diary off;
else
h = waitbar(0,'Taking subdraws...');
end
stock_param = zeros(MAX_nruns, npar);
stock_logpo = zeros(MAX_nruns,1);
@ -128,7 +126,7 @@ end
if options_.filter_step_ahead
stock_filter_step_ahead = zeros(naK,endo_nbr,gend+ ...
options_.filter_step_ahead(end),MAX_naK);
options_.filter_step_ahead(end),MAX_naK);
run_smoother = 1;
end
if options_.forecast
@ -139,153 +137,128 @@ end
%if moments_varendo
% stock_moments = cell(MAX_momentsno,1);
%end
for b=1:B
[deep, logpo] = GetOneDraw(type);
set_all_parameters(deep);
dr = resol(oo_.steady_state,0);
%if moments_varendo
% stock_moments{irun(8)} = compute_model_moments(dr,M_,options_);
%end
if run_smoother
[alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK] = ...
DsgeSmoother(deep,gend,Y,data_index,missing_value);
if options_.loglinear
stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ...
repmat(log(dr.ys(dr.order_var)),1,gend);
stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ...
repmat(log(dr.ys(dr.order_var)),1,gend);
else
stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ...
repmat(dr.ys(dr.order_var),1,gend);
stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ...
repmat(dr.ys(dr.order_var),1,gend);
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
if horizon
yyyy = alphahat(iendo,i_last_obs);
yf = forcst2a(yyyy,dr,zeros(horizon,exo_nbr));
if options_.prefilter == 1
yf(:,IdObs) = yf(:,IdObs)+repmat(bayestopt_.mean_varobs', ...
horizon+maxlag,1);
end
yf(:,IdObs) = yf(:,IdObs)+(gend+[1-maxlag:horizon]')*trend_coeff';
if options_.loglinear == 1
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]')* ...
trend_coeff',[1,1,1]);
if options_.loglinear == 1
yf1 = yf1 + repmat(log(SteadyState'),[horizon+maxlag,1,1]);
else
yf1 = yf1 + repmat(SteadyState',[horizon+maxlag,1,1]);
end
% Store the variable mandatory for local/remote parallel computing.
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';
localVars.typee=type;
localVars.run_smoother=run_smoother;
localVars.gend=gend;
localVars.Y=Y;
localVars.data_index=data_index;
localVars.missing_value=missing_value;
localVars.options_.varobs=options_.varobs;
localVars.irun=irun;
localVars.endo_nbr=endo_nbr;
localVars.nvn=nvn;
localVars.naK=naK;
localVars.horizon=horizon;
localVars.iendo=iendo;
if horizon
localVars.i_last_obs=i_last_obs;
localVars.IdObs=IdObs;
end
localVars.exo_nbr=exo_nbr;
localVars.maxlag=maxlag;
localVars.MAX_nsmoo=MAX_nsmoo;
localVars.MAX_ninno=MAX_ninno;
localVars.MAX_nerro = MAX_nerro;
if naK
localVars.MAX_naK=MAX_naK;
end
localVars.MAX_nruns=MAX_nruns;
if horizon
localVars.MAX_nforc1=MAX_nforc1;
localVars.MAX_nforc2=MAX_nforc2;
end
localVars.MAX_momentsno = MAX_momentsno;
localVars.ifil=ifil;
irun = irun + ones(7,1);
if irun(1) > MAX_nsmoo || b == B
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');
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');
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');
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');
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');
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');
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');
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;
% save([DirectoryName '/' M_.fname '_moments' int2str(ifil(8)) '.mat'],'stock');
% irun(8) = 1;
% end
if exist('OCTAVE_VERSION')
printf('Taking subdraws: %3.f%% done\r', b/B);
else
waitbar(b/B,h);
if strcmpi(type,'posterior'),
b=0;
while b<=B
b = b + 1;
[x(b,:), logpost(b,1)] = GetOneDraw(type);
end
end
if ~strcmpi(type,'prior'),
localVars.x=x;
localVars.logpost=logpost;
end
b=0;
if isnumeric(options_.parallel),% | isunix, % for the moment exclude unix platform from parallel implementation
[fout] = prior_posterior_statistics_core(localVars,1,B,0);
else
[nCPU, totCPU, nBlockPerCPU] = distributeJobs(options_.parallel, 1, B);
for j=1:totCPU-1,
nfiles = ceil(nBlockPerCPU(j)/MAX_nsmoo);
ifil(1,j+1) =ifil(1,j)+nfiles;
nfiles = ceil(nBlockPerCPU(j)/MAX_ninno);
ifil(2,j+1) =ifil(2,j)+nfiles;
nfiles = ceil(nBlockPerCPU(j)/MAX_nerro);
ifil(3,j+1) =ifil(3,j)+nfiles;
nfiles = ceil(nBlockPerCPU(j)/MAX_naK);
ifil(4,j+1) =ifil(4,j)+nfiles;
nfiles = ceil(nBlockPerCPU(j)/MAX_nruns);
ifil(5,j+1) =ifil(5,j)+nfiles;
nfiles = ceil(nBlockPerCPU(j)/MAX_nforc1);
ifil(6,j+1) =ifil(6,j)+nfiles;
nfiles = ceil(nBlockPerCPU(j)/MAX_nforc2);
ifil(7,j+1) =ifil(7,j)+nfiles;
% nfiles = ceil(nBlockPerCPU(j)/MAX_momentsno);
% ifil(8,j+1) =ifil(8,j)+nfiles;
end
localVars.ifil = ifil;
globalVars = struct('M_',M_, ...
'options_', options_, ...
'bayestopt_', bayestopt_, ...
'estim_params_', estim_params_, ...
'oo_', oo_);
% which files have to be copied to run remotely
NamFileInput(1,:) = {'',[M_.fname '_static.m']};
NamFileInput(2,:) = {'',[M_.fname '_dynamic.m']};
if options_.steadystate_flag,
NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '_steadystate.m']};
end
[fout] = masterParallel(options_.parallel, 1, B,NamFileInput,'prior_posterior_statistics_core', localVars,globalVars, options_.parallel_info);
end
ifil = fout(end).ifil;
if exist('OCTAVE_VERSION')
printf('\n');
diary on;
else
close(h)
if exist('h')
close(h)
end
end
stock_gend=gend;
stock_data=Y;
save([DirectoryName '/' M_.fname '_data.mat'],'stock_gend','stock_data');
if ~isnumeric(options_.parallel),
leaveSlaveOpen = options_.parallel_info.leaveSlaveOpen;
if options_.parallel_info.leaveSlaveOpen == 0,
options_.parallel_info.leaveSlaveOpen = 1; % force locally to leave open remote matlab sessions (repeated pm3 calls)
end
end
if options_.smoother
pm3(endo_nbr,gend,ifil(1),B,'Smoothed variables',...
'',varlist,'tit_tex',M_.endo_names,...
'',M_.endo_names(1:M_.orig_endo_nbr, :),'tit_tex',M_.endo_names,...
varlist,'SmoothedVariables',[M_.dname '/metropolis'],'_smooth');
pm3(exo_nbr,gend,ifil(2),B,'Smoothed shocks',...
'',M_.exo_names,'tit_tex',M_.exo_names,...
@ -315,4 +288,12 @@ if options_.forecast
pm3(endo_nbr,horizon+maxlag,ifil(6),B,'Forecasted variables (point)',...
'',varlist,'tit_tex',M_.endo_names,...
varlist,'PointForecast',[M_.dname '/metropolis'],'_forc_point');
end
end
if ~isnumeric(options_.parallel),
options_.parallel_info.leaveSlaveOpen = leaveSlaveOpen;
if leaveSlaveOpen == 0,
closeSlave(options_.parallel),
end
end

238
matlab/prior_posterior_statistics_core.m Normal file
View File

@ -0,0 +1,238 @@
function myoutput=prior_posterior_statistics_core(myinputs,fpar,B,whoiam, ThisMatlab)
if nargin<4,
whoiam=0;
end
struct2local(myinputs);
global options_ oo_ M_ bayestopt_ estim_params_
DirectoryName = CheckPath('metropolis');
RemoteFlag = 0;
if whoiam
ifil=ifil(:,whoiam);
waitbarString = ['Please wait... Bayesian (posterior) subdraws (' int2str(fpar) 'of' int2str(B) ')...'];
if Parallel(ThisMatlab).Local,
waitbarTitle=['Local '];
else
waitbarTitle=[Parallel(ThisMatlab).PcName];
RemoteFlag = 1;
end
fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo);
else
if exist('OCTAVE_VERSION')
diary off;
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 = {};
end
for b=fpar:B
% [deep, logpo] = GetOneDraw(typee);
% set_all_parameters(deep);
% dr = resol(oo_.steady_state,0);
if strcmpi(typee,'prior')
[deep, logpo] = GetOneDraw(typee);
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);
if options_.loglinear
stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ...
repmat(log(dr.ys(dr.order_var)),1,gend);
stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ...
repmat(log(dr.ys(dr.order_var)),1,gend);
else
stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ...
repmat(dr.ys(dr.order_var),1,gend);
stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ...
repmat(dr.ys(dr.order_var),1,gend);
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
if horizon
yyyy = alphahat(iendo,i_last_obs);
yf = forcst2a(yyyy,dr,zeros(horizon,exo_nbr));
if options_.prefilter == 1
yf(:,IdObs) = yf(:,IdObs)+repmat(bayestopt_.mean_varobs', ...
horizon+maxlag,1);
end
yf(:,IdObs) = yf(:,IdObs)+(gend+[1-maxlag:horizon]')*trend_coeff';
if options_.loglinear == 1
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]')* ...
trend_coeff',[1,1,1]);
if options_.loglinear == 1
yf1 = yf1 + repmat(log(SteadyState'),[horizon+maxlag,1,1]);
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);
%
% TempPath=DirectoryName;
% DirectoryNamePar='C:\dynare_calcs\ModelTest\ls2003\metropolis'
% DirectoryName=DirectoryNamePar;
if irun(1) > MAX_nsmoo || b == B
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,
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']}];
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']}];
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']}];
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']}];
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']}];
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']}];
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;
% save([DirectoryName '/' M_.fname '_moments' int2str(ifil(8)) '.mat'],'stock');
% if RemoteFlag==1,
% OutputFileName_moments = [OutputFileName_moments; {[DirectoryName filesep], [M_.fname '_moments' int2str(ifil(8)) '.mat']}];
% end
% irun(8) = 1;
% end
% DirectoryName=TempPath;
if exist('OCTAVE_VERSION')
printf('Taking subdraws: %3.f%% done\r', b/B*100);
elseif ~whoiam,
waitbar(b/B,h);
end
if mod(b,10)==0 & whoiam,
fprintf('Done! \n');
waitbarString = [ 'Subdraw ' int2str(b) '/' int2str(B) ' done.'];
fMessageStatus((b-fpar+1)/(B-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo)
end
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];
end

4
matlab/random_walk_metropolis_hastings.m
View File

@ -69,7 +69,7 @@ localVars.varargin=varargin;
% tic,
if isnumeric(options_.parallel),
if isnumeric(options_.parallel) || (nblck-fblck)==0,
fout = random_walk_metropolis_hastings_core(localVars, fblck, nblck, 0);
record = fout.record;
@ -94,7 +94,7 @@ else
% from where to get back results
% NamFileOutput(1,:) = {[M_.dname,'/metropolis/'],'*.*'};
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'random_walk_metropolis_hastings_core', localVars, globalVars);
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'random_walk_metropolis_hastings_core', localVars, globalVars, options_.parallel_info);
for j=1:totCPU,
offset = sum(nBlockPerCPU(1:j-1))+fblck-1;
record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)));