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

TeX formatting and various things... 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@656 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
adjemian 2006-03-06 08:38:11 +00:00
parent 00fac4544c
commit 58cb6e42b6
3 changed files with 32 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
matlab/dynare_estimation.m
View File

@ -191,8 +191,7 @@ rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet
k = []; k = [];
k1 = []; k1 = [];
for i=1:n_varobs for i=1:n_varobs
k = [k strmatch(deblank(options_.varobs(i,:)),M_.endo_names(dr.order_var,:), ... k = [k strmatch(deblank(options_.varobs(i,:)),M_.endo_names(dr.order_var,:),'exact')];
'exact')];
k1 = [k1 strmatch(deblank(options_.varobs(i,:)),M_.endo_names, 'exact')]; k1 = [k1 strmatch(deblank(options_.varobs(i,:)),M_.endo_names, 'exact')];
end end
@ -257,29 +256,9 @@ if options_.mode_compute > 0 & options_.posterior_mode_estimation
disp(sprintf('Objective function at mode: %f',fval)) disp(sprintf('Objective function at mode: %f',fval))
disp(sprintf('Objective function at mode: %f',DsgeLikelihood(xparam1,gend,data))) disp(sprintf('Objective function at mode: %f',DsgeLikelihood(xparam1,gend,data)))
elseif options_.mode_compute == 5 elseif options_.mode_compute == 5
if isfield(options_,'hess') flag = 0;
flag = options_.hess; [xparam1, hh, gg, fval] = newrat('DsgeLikelihood',xparam1,[],[],flag,gend,data);
else eval(['save ' M_.fname '_mode xparam1 hh gg fval;']);
flag = 1;
end
if ~exist('igg'), % by M. Ratto
hh=[];
gg=[];
igg=[];
end % by M. Ratto
if isfield(options_,'ftol')
crit = options_.ftol;
else
crit = 1.e-7;
end
if isfield(options_,'nit')
nit = options_.nit;
else
nit=1000;
end
%[xparam1, hh, gg, fval] = newrat('DsgeLikelihood',xparam1,[],[],flag,gend,data);
[xparam1, hh, gg, fval, invhess] = newrat('DsgeLikelihood',xparam1,hh,gg,igg,crit,nit,flag,gend,data);
eval(['save ' M_.fname '_mode xparam1 hh gg fval invhess;']);
end end
if options_.mode_compute ~= 5 if options_.mode_compute ~= 5
hh = reshape(hessian('DsgeLikelihood',xparam1,gend,data),nx,nx); hh = reshape(hessian('DsgeLikelihood',xparam1,gend,data),nx,nx);
@ -483,9 +462,11 @@ elseif ~any(bayestopt_.pshape > 0) & options_.posterior_mode_estimation
end end
OutputDirectoryName = CheckPath('Output');
if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior mode) Latex output if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior mode) Latex output
if np if np
filename = [M_.fname '_Posterior_Mode_1.TeX']; filename = [OutputDirectoryName '\' M_.fname '_Posterior_Mode_1.TeX'];
fidTeX = fopen(filename,'w'); fidTeX = fopen(filename,'w');
fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n'); fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (parameters)\n'); fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (parameters)\n');
@ -502,7 +483,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
ip = nvx+nvn+ncx+ncn+1; ip = nvx+nvn+ncx+ncn+1;
for i=1:np for i=1:np
fprintf(fidTeX,'$%s$ & %s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n',... fprintf(fidTeX,'$%s$ & %s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n',...
M_.param_names_tex(estim_params_.param_vals(i,1),:),...%deblank(estim_params_.tex(i,:)),... M_.param_names_tex(estim_params_.param_vals(i,1),:),...
deblank(pnames(bayestopt_.pshape(ip)+1,:)),... deblank(pnames(bayestopt_.pshape(ip)+1,:)),...
bayestopt_.pmean(ip),... bayestopt_.pmean(ip),...
estim_params_.param_vals(i,6),... estim_params_.param_vals(i,6),...
@ -520,7 +501,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fclose(fidTeX); fclose(fidTeX);
end end
if nvx if nvx
TeXfile = [M_.fname '_Posterior_Mode_2.TeX']; TeXfile = [OutputDirectoryName '\' M_.fname '_Posterior_Mode_2.TeX'];
fidTeX = fopen(TeXfile,'w'); fidTeX = fopen(TeXfile,'w');
fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n'); fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (standard deviation of structural shocks)\n'); fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (standard deviation of structural shocks)\n');
@ -556,7 +537,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fclose(fidTeX); fclose(fidTeX);
end end
if nvn if nvn
TeXfile = [M_.fname '_Posterior_Mode_3.TeX']; TeXfile = [OutputDirectoryName '\' M_.fname '_Posterior_Mode_3.TeX'];
fidTeX = fopen(TeXfile,'w'); fidTeX = fopen(TeXfile,'w');
fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n'); fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (standard deviation of measurement errors)\n'); fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (standard deviation of measurement errors)\n');
@ -590,7 +571,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fclose(fidTeX); fclose(fidTeX);
end end
if ncx if ncx
TeXfile = [M_.fname '_Posterior_Mode_4.TeX']; TeXfile = [OutputDirectoryName '\' M_.fname '_Posterior_Mode_4.TeX'];
fidTeX = fopen(TeXfile,'w'); fidTeX = fopen(TeXfile,'w');
fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n'); fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (correlation of structural shocks)\n'); fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (correlation of structural shocks)\n');
@ -625,7 +606,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fclose(fidTeX); fclose(fidTeX);
end end
if ncn if ncn
TeXfile = [M_.fname '_Posterior_Mode_5.TeX']; TeXfile = [OutputDirectoryName '\' M_.fname '_Posterior_Mode_5.TeX'];
fidTeX = fopen(TeXfile,'w'); fidTeX = fopen(TeXfile,'w');
fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n'); fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (correlation of measurement errors)\n'); fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (correlation of measurement errors)\n');

41
matlab/metropolis.m
View File

@ -4,7 +4,7 @@ global M_ oo_ options_ bayestopt_ estim_params_
bayestopt_.penalty = 1e8; bayestopt_.penalty = 1e8;
DirectoryName = CheckPath('metropolis'); MhDirectoryName = CheckPath('metropolis');
nblck = options_.mh_nblck; nblck = options_.mh_nblck;
nruns = ones(nblck,1)*options_.mh_replic; nruns = ones(nblck,1)*options_.mh_replic;
@ -23,9 +23,9 @@ if options_.load_mh_file == 0
disp('MH: One Chain mode.') disp('MH: One Chain mode.')
end end
% Delete old mh files... % Delete old mh files...
files = eval(['dir(''' DirectoryName '/' M_.fname '_mh*_blck*.mat'');']); files = dir([ MhDirectoryName '/' M_.fname '_mh*_blck*.mat']);
if length(files) if length(files)
delete([ DirectoryName '/' M_.fname '_mh*_blck*.mat']); delete([ MhDirectoryName '/' M_.fname '_mh*_blck*.mat']);
disp('MH: Old _mh files succesfully erased!') disp('MH: Old _mh files succesfully erased!')
end end
% Initial values... % Initial values...
@ -78,9 +78,9 @@ if options_.load_mh_file == 0
fline = ones(nblck,1); fline = ones(nblck,1);
NewFile = ones(nblck,1); NewFile = ones(nblck,1);
% Creation of the mh-history file: % Creation of the mh-history file:
file = eval(['dir(''' DirectoryName '/' M_.fname '_mh_history.mat'');']); file = dir([MhDirectoryName '/' M_.fname '_mh_history.mat']);
if length(files) if length(files)
delete([ DirectoryName '/' M_.fname '_mh_history.mat']); delete([ MhDirectoryName '/' M_.fname '_mh_history.mat']);
disp('MH: Old mh_history file succesfully erased!') disp('MH: Old mh_history file succesfully erased!')
end end
AnticipatedNumberOfFiles = floor(nruns(1)/MAX_nruns); AnticipatedNumberOfFiles = floor(nruns(1)/MAX_nruns);
@ -102,8 +102,8 @@ if options_.load_mh_file == 0
save([DirectoryName '/' M_.fname '_mh_history'],'record'); save([DirectoryName '/' M_.fname '_mh_history'],'record');
elseif options_.load_mh_file == 1% Here we consider previous mh files (previous mh did not crash). elseif options_.load_mh_file == 1% Here we consider previous mh files (previous mh did not crash).
disp('MH: I''m loading past metropolis-hastings simulations...') disp('MH: I''m loading past metropolis-hastings simulations...')
file = eval(['dir(''' DirectoryName '/' M_.fname '_mh_history.mat'');']); file = dir([ MhDirectoryName '/' M_.fname '_mh_history.mat' ]);
files = eval(['dir(''' DirectoryName '/' M_.fname '_mh*.mat'');']); files = dir([ MhDirectoryName '/' M_.fname '_mh*.mat']);
if ~length(files) if ~length(files)
disp('MH:: FAILURE! there is no MH file to load here!') disp('MH:: FAILURE! there is no MH file to load here!')
return return
@ -112,7 +112,7 @@ elseif options_.load_mh_file == 1% Here we consider previous mh files (previous
disp('MH:: FAILURE! there is no MH-history file!') disp('MH:: FAILURE! there is no MH-history file!')
return return
else else
load([ DirectoryName '/' M_.fname '_mh_history']) load([ MhDirectoryName '/' M_.fname '_mh_history'])
end end
past_number_of_blocks = record.Nblck; past_number_of_blocks = record.Nblck;
if past_number_of_blocks ~= nblck if past_number_of_blocks ~= nblck
@ -148,7 +148,7 @@ elseif options_.load_mh_file == 1% Here we consider previous mh files (previous
record.MhDraws(end,3) = AnticipatedNumberOfLinesInTheLastFile; record.MhDraws(end,3) = AnticipatedNumberOfLinesInTheLastFile;
randn('state',record.Seeds.Normal); randn('state',record.Seeds.Normal);
rand('state',record.Seeds.Unifor); rand('state',record.Seeds.Unifor);
save([DirectoryName '/' M_.fname '_mh_history'],'record'); save([MhDirectoryName '/' M_.fname '_mh_history'],'record');
disp(['MH: ... It''s done. I''ve loaded ' int2str(NumberOfPreviousSimulations) ' simulations.']) disp(['MH: ... It''s done. I''ve loaded ' int2str(NumberOfPreviousSimulations) ' simulations.'])
disp(' ') disp(' ')
elseif options_.load_mh_file == -1% The previous metropolis-hastings elseif options_.load_mh_file == -1% The previous metropolis-hastings
@ -156,13 +156,13 @@ elseif options_.load_mh_file == -1% The previous metropolis-hastings
% recover the saved draws... % recover the saved draws...
disp('MH: Recover mode!') disp('MH: Recover mode!')
disp(' ') disp(' ')
file = eval(['dir(''' DirectoryName '/' M_.fname '_mh_history.mat'');']); file = dir([MhDirectoryName '/' M_.fname '_mh_history.mat']);
if ~length(file) if ~length(file)
disp('MH:: FAILURE! there is no MH-history file!') disp('MH:: FAILURE! there is no MH-history file!')
return return
else else
load([ DirectoryName '/' M_.fname '_mh_history']) load([ MhDirectoryName '/' M_.fname '_mh_history'])
end end
nblck = record.Nblck; nblck = record.Nblck;
options_.mh_nblck = nblck; options_.mh_nblck = nblck;
if size(record.MhDraws,1) == 1 if size(record.MhDraws,1) == 1
@ -202,12 +202,12 @@ elseif options_.load_mh_file == -1% The previous metropolis-hastings
ExpectedNumberOfMhFilesPerBlock = sum(record.MhDraws(:,2),1); ExpectedNumberOfMhFilesPerBlock = sum(record.MhDraws(:,2),1);
ExpectedNumberOfMhFiles = ExpectedNumberOfMhFilesPerBlock*nblck; ExpectedNumberOfMhFiles = ExpectedNumberOfMhFilesPerBlock*nblck;
% I count the total number of saved mh files... % I count the total number of saved mh files...
AllMhFiles = eval(['dir(''' DirectoryName '/' M_.fname '_mh*_blck*.mat'');']); AllMhFiles = dir([MhDirectoryName '/' M_.fname '_mh*_blck*.mat']);
TotalNumberOfMhFiles = length(AllMhFiles); TotalNumberOfMhFiles = length(AllMhFiles);
% I count the number of saved mh files per block % I count the number of saved mh files per block
NumberOfMhFilesPerBlock = zeros(nblck,1); NumberOfMhFilesPerBlock = zeros(nblck,1);
for i = 1:nblck for i = 1:nblck
BlckMhFiles = eval(['dir(''' DirectoryName '/' M_.fname '_mh*_blck' int2str(i) '.mat'');']); BlckMhFiles = dir([ MhDirectoryName '/' M_.fname '_mh*_blck' int2str(i) '.mat']);
NumberOfMhFilesPerBlock(i) = length(BlckMhFiles); NumberOfMhFilesPerBlock(i) = length(BlckMhFiles);
end end
tmp = NumberOfMhFilesPerBlock(1); b = 1; tmp = NumberOfMhFilesPerBlock(1); b = 1;
@ -229,8 +229,7 @@ elseif options_.load_mh_file == -1% The previous metropolis-hastings
% (if there was a complete session before the crash) ? % (if there was a complete session before the crash) ?
if OldMh if OldMh
ante = sum(record.MhDraws(1:end-1,2),1); ante = sum(record.MhDraws(1:end-1,2),1);
load(['./' DirectoryName '/' M_.fname '_mh' int2str(ante) '_blck' ... load([MhDirectoryName '/' M_.fname '_mh' int2str(ante) '_blck' int2str(CrashedBlck) '.mat'],'logpo2');
int2str(CrashedBlck) '.mat'],'logpo2');
if length(logpo2) == MAX_nruns if length(logpo2) == MAX_nruns
IsTheLastFileOfThePreviousMhFull = 1; IsTheLastFileOfThePreviousMhFull = 1;
else else
@ -245,7 +244,7 @@ elseif options_.load_mh_file == -1% The previous metropolis-hastings
MhFileNumber = ante; MhFileNumber = ante;
while MhFileExist while MhFileExist
MhFileNumber = MhFileNumber + 1; MhFileNumber = MhFileNumber + 1;
if ~exist(['./' DirectoryName '/' M_.fname '_mh' int2str(MhFileNumber) '_blck' int2str(CrashedBlck) '.mat']) if ~exist([MhDirectoryName '/' M_.fname '_mh' int2str(MhFileNumber) '_blck' int2str(CrashedBlck) '.mat'])
MhFileExist = 0; MhFileExist = 0;
end end
end end
@ -267,7 +266,7 @@ elseif options_.load_mh_file == -1% The previous metropolis-hastings
nruns(CrashedBlck) = nruns(CrashedBlck)-NumberOfSavedDraws; nruns(CrashedBlck) = nruns(CrashedBlck)-NumberOfSavedDraws;
% I initialize with the last saved mh file of the inccomplete % I initialize with the last saved mh file of the inccomplete
% block: % block:
load(['./' DirectoryName '/' M_.fname '_mh' int2str(MhFileNumber-1) '_blck' int2str(CrashedBlck) '.mat']); load([MhDirectoryName '/' M_.fname '_mh' int2str(MhFileNumber-1) '_blck' int2str(CrashedBlck) '.mat']);
ilogpo2(CrashedBlck) = logpo2(end); ilogpo2(CrashedBlck) = logpo2(end);
ix2(CrashedBlck,:) = x2(end,:); ix2(CrashedBlck,:) = x2(end,:);
NewFile(CrashedBlck) = MhFileNumber; NewFile(CrashedBlck) = MhFileNumber;
@ -319,7 +318,7 @@ for b = fblck:nblck
prtfrc = j/nruns(b); prtfrc = j/nruns(b);
waitbar(prtfrc,hh,[ '(' int2str(b) '/' int2str(nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]); waitbar(prtfrc,hh,[ '(' int2str(b) '/' int2str(nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]);
if (irun == InitSizeArray(b)) | (j == nruns(b)) % Now I save the simulations if (irun == InitSizeArray(b)) | (j == nruns(b)) % Now I save the simulations
save([DirectoryName '/' M_.fname '_mh' int2str(NewFile(b)) '_blck' int2str(b)],'x2','logpo2'); save([MhDirectoryName '/' M_.fname '_mh' int2str(NewFile(b)) '_blck' int2str(b)],'x2','logpo2');
InitSizeArray(b) = min(nruns(b)-j,MAX_nruns); InitSizeArray(b) = min(nruns(b)-j,MAX_nruns);
if j == nruns(b) % I record the last draw... if j == nruns(b) % I record the last draw...
record.LastParameters(b,:) = x2(end,:); record.LastParameters(b,:) = x2(end,:);
@ -342,8 +341,8 @@ for b = fblck:nblck
end% End of the loop over the mh-blocks. end% End of the loop over the mh-blocks.
record.Seeds.Normal = randn('state'); record.Seeds.Normal = randn('state');
record.Seeds.Unifor = rand('state'); record.Seeds.Unifor = rand('state');
save([DirectoryName '/' M_.fname '_mh_history'],'record'); save([MhDirectoryName '/' M_.fname '_mh_history'],'record');
disp(['MH: Number of mh files : ' int2str(NewFile(1)) ' per block.']) disp(['MH: Number of mh files : ' int2str(NewFile(1)) ' per block.'])
disp(['MH: Total number of generated files : ' int2str(NewFile(1)*nblck) '.']) disp(['MH: Total number of generated files : ' int2str(NewFile(1)*nblck) '.'])
disp(['MH: Total number of iterations : ' int2str((NewFile(1)-1)*MAX_nruns+irun-1) '.']) disp(['MH: Total number of iterations : ' int2str((NewFile(1)-1)*MAX_nruns+irun-1) '.'])
disp(' ') disp(' ')

1
matlab/posterior_moments.m
View File

@ -39,6 +39,5 @@ else
bandwidth = 0; % Rule of thumb optimal bandwidth parameter. bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton. kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton.
optimal_bandwidth = mh_optimal_bandwidth(xx,length(xx),bandwidth,kernel_function); optimal_bandwidth = mh_optimal_bandwidth(xx,length(xx),bandwidth,kernel_function);
% [abscissa,f] = kernel_density_estimate(,,,)
[density(:,1),density(:,2)] = kernel_density_estimate(xx,number_of_grid_points,optimal_bandwidth,kernel_function); [density(:,1),density(:,2)] = kernel_density_estimate(xx,number_of_grid_points,optimal_bandwidth,kernel_function);
end end