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

re-organisation of folders to package separately GSA tools 

git-svn-id: file:///home/sebastien/dynare/gsa_dyn@5 f1850c17-3b45-254b-b221-fcb05880fee1
time-shift
rattoma 2006-11-22 08:20:49 +00:00
parent 5fdf16bdcc
commit 554c3c09ed
65 changed files with 0 additions and 7280 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

118
matlab/LoadPosteriorSmoother.m
View File

@ -1,118 +0,0 @@
function PosteriorSmoother(type)
% stephane.adjemian@ens.fr [09-25-2005]
global options_ estim_params_ oo_ M_
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
offset = npar-np;
naK = length(options_.filter_step_ahead);
%%
endo_nbr=M_.endo_nbr;
exo_nbr=M_.exo_nbr;
nvobs = size(options_.varobs,1);
%%
CheckPath('Plots/');
DirectoryName = CheckPath('metropolis');
%%
%%
varlist = options_.varlist;
if isempty(varlist)
varlist = M_.endo_names;
SelecVariables = transpose(1:M_.endo_nbr);
nvar = M_.endo_nbr;
else
nvar = size(varlist,1);
SelecVariables = [];
for i=1:nvar
if ~isempty(strmatch(varlist(i,:),M_.endo_names,'exact'))
SelecVariables = [SelecVariables;strmatch(varlist(i,:),M_.endo_names,'exact')];
end
end
end
h = waitbar(0,'Load Bayesian smoother...');
filsmooth = dir([DirectoryName '/' M_.fname '_smooth*.mat']);
y0=[];
for j=1:length(filsmooth),
load([DirectoryName '/' M_.fname '_smooth',num2str(j),'.mat']);
if isempty(y0),
y0=stock;
else
%y0(:,:,size(y0,3):size(y0,3)+size(stock,3))=stock;
y0=cat(3,y0,stock);
end
end
if nvx,
filinno = dir([DirectoryName '/' M_.fname '_inno*.mat']);
end
if nvn
filerror = dir([DirectoryName '/' M_.fname '_error*.mat']);
end
if naK
filfilt = dir([DirectoryName '/' M_.fname '_filter*.mat']);
end
filparam = dir([DirectoryName '/' M_.fname '_param*.mat']);
for b=1:B
%deep = GetOneDraw(NumberOfDraws,FirstMhFile,LastMhFile,FirstLine,MAX_nruns,DirectoryName);
deep = GetOneDraw(type);
if irun1 > MAX_nsmoo | b == B
if b == B
stock_smooth = stock_smooth(:,:,1:irun1-1);
end
stock = stock_smooth;
save([DirectoryName '/' M_.fname '_smooth' int2str(ifil1)],'stock');
ifil1 = ifil1 + 1;
irun1 = 1;
end
if nvx & (irun2 > MAX_ninno | b == B)
if b == B
stock_innov = stock_innov(:,:,1:irun2-1);
end
stock = stock_innov;
save([DirectoryName '/' M_.fname '_inno' int2str(ifil2)],'stock');
ifil2 = ifil2 + 1;
irun2 = 1;
end
if nvn & (irun3 > MAX_error | b == B)
if b == B
stock_error = stock_error(:,:,1:irun3-1);
end
stock = stock_error;
save([DirectoryName '/' M_.fname '_error' int2str(ifil3)],'stock');
ifil3 = ifil3 + 1;
irun3 = 1;
end
if naK & (irun4 > MAX_naK | b == B)
if b == B
stock_filter = stock_filter(:,:,:,1:irun4-1);
end
stock = stock_filter;
save([DirectoryName '/' M_.fname '_filter' int2str(ifil4)],'stock');
ifil4 = ifil4 + 1;
irun4 = 1;
end
if irun5 > MAX_nruns | b == B
if b == B
stock_param = stock_param(1:irun5-1,:);
end
stock = stock_param;
save([DirectoryName '/' M_.fname '_param' int2str(ifil5)],'stock');
ifil5 = ifil5 + 1;
irun5 = 1;
end
waitbar(b/B,h);
end
close(h)

100
matlab/PlotPosteriorHistograms.m
View File

@ -1,100 +0,0 @@
function PlotPosteriorHistograms()
% stephane.adjemian@ens.fr [09-09-2005]
global estim_params_ M_ options_ bayestopt_ oo_
OutputDirectoryName = CheckPath('Output');
TeX = options_.TeX;
nblck = options_.mh_nblck;
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
MaxNumberOfPlotPerFigure = 9;% The square root must be an integer!
nn = sqrt(MaxNumberOfPlotPerFigure);
figurename = 'Priors and posteriors hitograms';
if TeX
fidTeX = fopen([OutputDirectoryName '/' M_.fname '_PriorsAndPosteriorsHist.TeX'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by PlotPosteriorHistograms.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
end
figunumber = 0;
subplotnum = 0;
for i=1:npar
subplotnum = subplotnum+1;
if subplotnum == 1
figunumber = figunumber+1;
if options_.nograph
hfig = figure('Name',figurename,'Visible','off');
else
hfig = figure('Name',figurename);
end
end
if subplotnum == 1
if TeX
TeXNAMES = [];
end
NAMES = [];
end
[nam,texnam] = get_the_name(i,TeX);
NAMES = strvcat(NAMES,nam);
if TeX
TeXNAMES = strvcat(TeXNAMES,texnam);
end
hist
%[x2,f2,abscissa,dens,binf2,bsup2] = draw_prior_density(i);
top2 = max(f2);
top1 = max(f1);
top0 = max([top1;top2]);
binf1 = x1(1);
bsup1 = x1(end);
borneinf = min(binf1,binf2);
bornesup = max(bsup1,bsup2);
subplot(nn,nn,subplotnum)
hh = plot(x2,f2,'-k','linewidth',2);
set(hh,'color',[0.7 0.7 0.7]);
hold on;
plot(x1,f1,'-k','linewidth',2);
plot( [pmode pmode], [0.0 1.1*top0], '--g', 'linewidth', 2);
box on;
axis([borneinf bornesup 0 1.1*top0]);
title(nam,'Interpreter','none');
hold off;
drawnow
if subplotnum == MaxNumberOfPlotPerFigure | i == npar;
eval(['print -depsc2 ' OutputDirectoryName '/' M_.fname '_PriorsAndPosteriors' int2str(figunumber)]);
eval(['print -dpdf ' OutputDirectoryName '/' M_.fname '_PriorsAndPosteriors' int2str(figunumber)]);
if options_.nograph,
set(hfig,'Visible','on');
end
saveas(hfig,[OutputDirectoryName '/' M_.fname '_PriorsAndPosteriors' int2str(figunumber) '.fig']);
if TeX
fprintf(fidTeX,'\\begin{figure}[H]\n');
for j = 1:size(NAMES,1)
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(j,:)),deblank(TeXNAMES(j,:)));
end
fprintf(fidTeX,'\\centering\n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_PriorsAndPosteriors%s}\n',M_.fname,int2str(figunumber));
fprintf(fidTeX,'\\caption{Priors and posteriors.}');
fprintf(fidTeX,'\\label{Fig:PriorsAndPosteriors:%s}\n',int2str(figunumber));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
if i == npar
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
end
end
if options_.nograph,
close(hfig),
end
subplotnum = 0;
end
end

143
matlab/PosteriorParametersHistograms.m
View File

@ -1,143 +0,0 @@
function PosteriorParametersHistograms()
% stephane.adjemian@ens.fr [09-09-2005]
global estim_params_ M_ options_ bayestopt_ oo_
TeX = options_.TeX;
nblck = options_.mh_nblck;
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
nx = nvx+nvn+ncx+ncn+np;
DirectoryName = CheckPath('metropolis');
OutputDirectoryName = CheckPath('Output');
load([ DirectoryName '/' M_.fname '_mh_history'])
FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2))
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
FirstMhFile = record.KeepedDraws.FirstMhFile;
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
clear record;
MaxNumberOfPlotPerFigure = 9;% The square root must be an integer!
nn = sqrt(MaxNumberOfPlotPerFigure);
figurename = 'Priors and posteriors hitograms';
if TeX
fidTeX = fopen([OutputDirectoryName '/' M_.fname '_PosteriorsHist.TeX'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by PlotPosteriorHistograms.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
end
figunumber = 0;
subplotnum = 0;
for ip = 1:nx,
subplotnum = subplotnum+1;
if subplotnum == 1
figunumber = figunumber+1;
if options_.nograph
hfig = figure('Name',figurename,'Visible','off');
else
hfig = figure('Name',figurename);
end
end
if subplotnum == 1
if TeX
TeXNAMES = [];
end
NAMES = [];
end
[nam,texnam] = get_the_name(ip,TeX);
NAMES = strvcat(NAMES,nam);
if TeX
TeXNAMES = strvcat(TeXNAMES,texnam);
end
if ip<=nvx
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
%[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
[N,X] = hist(Draws,30);
k = estim_params_.var_exo(ip,1);
name = deblank(M_.exo_names(k,:));
eval(['pmode = oo_.posterior_mode.shocks_std.' name ';'])
elseif ip<=nvn+nvx
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
%[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
[N,X] = hist(Draws,30);
name = deblank(options_.varobs(estim_params_.var_endo(ip-nvx,1),:));
eval(['pmode = oo_.posterior_mode.measurement_errors_std.' name ';'])
elseif ip<=nvn+nvx+ncx
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
%[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
[N,X] = hist(Draws,30);
k1 = estim_params_.corrx(ip-nvx-nvn,1);
k2 = estim_params_.corrx(ip-nvx-nvn,2);
name = [deblank(M_.exo_names(k1,:)) ',' deblank(M_.exo_names(k2,:))];
NAME = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
eval(['pmode = oo_.posterior_mode.shocks_corr.' name ';'])
elseif ip<=nvn+nvx+ncx+ncn
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
%[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
[N,X] = hist(Draws,30);
k1 = estim_params_.corrn(ip-(nvx+nvn+ncx),1);
k2 = estim_params_.corrn(ip-(nvx+nvn+ncx),2);
name = [deblank(M_.endo_names(k1,:)) ',' deblank(M_.endo_names(k2,:))];
NAME = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
eval(['pmode = oo_.posterior_mode.measurement_errors_corr.' name ';'])
else
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
%[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
[N,X] = hist(Draws,30);
name = bayestopt_.name{ip};
eval(['pmode = oo_.posterior_mode.parameters.' name ';'])
end
subplot(nn,nn,subplotnum)
f=N./length(Draws);
top0=max(f);
hh = bar(X,f);
minbound=min(Draws);
maxbound=max(Draws);
hold on;
plot( [pmode pmode], [0.0 1.1*top0], '--g', 'linewidth', 2);
box on;
axis([minbound maxbound 0 1.1*top0]);
title(nam,'Interpreter','none');
hold off;
drawnow
if subplotnum == MaxNumberOfPlotPerFigure | ip == nx;
eval(['print -depsc2 ' OutputDirectoryName '/' M_.fname '_PosteriorHist' int2str(figunumber)]);
eval(['print -dpdf ' OutputDirectoryName '/' M_.fname '_PosteriorHist' int2str(figunumber)]);
if options_.nograph,
set(hfig,'Visible','on');
end
saveas(hfig,[OutputDirectoryName '/' M_.fname '_PosteriorHist' int2str(figunumber) '.fig']);
if TeX
fprintf(fidTeX,'\\begin{figure}[H]\n');
for j = 1:size(NAMES,1)
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(j,:)),deblank(TeXNAMES(j,:)));
end
fprintf(fidTeX,'\\centering\n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_PosteriorHist%s}\n',M_.fname,int2str(figunumber));
fprintf(fidTeX,'\\caption{Priors and posteriors.}');
fprintf(fidTeX,'\\label{Fig:PosteriorHist:%s}\n',int2str(figunumber));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
if i == npar
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
end
end
if options_.nograph,
close(hfig),
end
subplotnum = 0;
end
end

10
matlab/clear_dynare_project.m
View File

@ -1,10 +0,0 @@
dname=ls('*.mod');
dname=dname(1:end-4);
delete('*.eps')
delete('*.pdf')
delete('*.asv')
delete('*.fig')
delete('*.tex')
delete('*.log*')
flag=rmdir(dname,'s');

59
matlab/cumdens.m
View File

@ -1,59 +0,0 @@
function [xcum, paracdf] = cumdens(para, pshape, p1, p2, p3, p4)
% This procedure transforms x vectors into cumulative values
% pshape: 0 is point mass, both para and p2 are ignored
% 1 is BETA(mean,stdd)
% 2 is GAMMA(mean,stdd)
% 3 is NORMAL(mean,stdd)
% 4 is INVGAMMA(s^2,nu)
% 5 is UNIFORM [p1,p2]
% Adapted by M. Ratto from MJ priordens.m
nprio = length(pshape);
i = 1;
while i <= nprio;
a = 0;
b = 0;
if pshape(i) == 1; % (generalized) BETA Prior
mu = (p1(i)-p3(i))/(p4(i)-p3(i));
stdd = p2(i)/(p4(i)-p3(i));
a = (1-mu)*mu^2/stdd^2 - mu;
b = a*(1/mu - 1);
%lnprior = lnprior + lpdfgbeta(para(i),a,b,p3(i),p4(i)) ;
para(:,i) = (para(:,i)-p3(i))./(p4(i)-p3(i));
xcum(:,i) = betacdf(para(:,i),a,b) ;
as=num2str(a,'%20.16e');
bs=num2str(b,'%20.16e');
p3s=num2str(p3(i),'%20.16e');
p4s=num2str(p4(i),'%20.16e');
paracdf{i} = inline(['betacdf((X-',p3s,')./(',p4s,'-',p3s,'),',as,',',bs,')'],'X');
elseif pshape(i) == 2; % GAMMA PRIOR
b = p2(i)^2/(p1(i)-p3(i));
a = (p1(i)-p3(i))/b;
%lnprior = lnprior + lpdfgam(para(i)-p3(i),a,b);
xcum(:,i) = gamcdf(para(:,i)-p3(i),a,b);
as=num2str(a,'%20.16e');
bs=num2str(b,'%20.16e');
p3s=num2str(p3(i),'%20.16e');
paracdf{i} = inline(['gamcdf((X-',p3s,'),',as,',',bs,')'],'X');
elseif pshape(i) == 3; % GAUSSIAN PRIOR
%lnprior = lnprior + lpdfnorm(para(i),p1(i),p2(i));
xcum(:,i) = normcdf(para(:,i),p1(i),p2(i));
p1s=num2str(p1(i),'%20.16e');
p2s=num2str(p2(i),'%20.16e');
paracdf{i} = inline(['normcdf(X,',p1s,',',p2s,')'],'X');
elseif pshape(i) == 4; % INVGAMMA1 PRIOR
%lnprior = lnprior + lpdfig1(para(i),p1(i),p2(i));
xcum(:,i)=para(:,i);
elseif pshape(i) == 5; % UNIFORM PRIOR
%lnprior = lnprior + log(1/(p2(i)-p1(i)));
xcum(:,i) = (para(:,i)-p1(i))./(p2(i)-p1(i));
p1s=num2str(p1(i),'%20.16e');
p2s=num2str(p2(i),'%20.16e');
paracdf{i} = inline(['(X-',p1s,')./(',p2s,'-',p1s,')'],'X');
% elseif pshape(i) == 6; % INVGAMMA2 PRIOR
% lnprior = lnprior + lpdfig2(para(i),p1(i),p2(i));
end;
i = i+1;
end;

95
matlab/dirf_comp.m
View File

@ -1,95 +0,0 @@
function dirf_comp(var_list_, s1, s2, vargin, texname)
global M_ options_
DirectoryName = CheckPath('Output');
fname_ = M_.fname;
if nargin<5,
for j=1:M_.endo_nbr,
texname{j}=deblank(M_.endo_names_tex(j,:));
end
end
iendo=[];
for i=1:size(var_list_,1)
iendo=[iendo, strmatch(deblank(var_list_(i,:)),M_.endo_names,'exact')];
if nargin<5,
if strmatch('\log',texname{iendo(i)});
texname{iendo(i)}=texname{iendo(i)}(2:end);
end
else
if strmatch('\log',texname{i});
texname{i}=texname{i}(2:end);
end
end
end
for i=1:length(vargin)
nfigsav=0;
nfig=0;
nplo=0;
for j=1:size(var_list_,1),
if isfield(s1,[deblank(var_list_(j,:)),'_',vargin{i}]),
nplo=nplo+1;
if mod(nplo,9)==1,
figure('name',['Comparison of orthogonalised shocks to ',vargin{i}])
nfig=nfig+1;
end
subplot(3,3,nplo-9*(nfig-1))
y=getfield(s1,[deblank(var_list_(j,:)),'_',vargin{i}]);
if isfield(s2,[deblank(var_list_(j,:)),'_',vargin{i}])
y=[y getfield(s2,[deblank(var_list_(j,:)),'_',vargin{i}])];
else
y=[y zeros(size(y))];
end
if ~exist('ys');
ys=zeros([size(y), size(var_list_,1)]);
end
ys(:,:,j)=ys(:,:,j)+y;
x=[1:length(y)];
plot(x,y(:,1),'k',x,y(:,2),':k'),
if options_.TeX,
if nargin<5
title(deblank(texname{iendo(j)}),'interpreter','tex')
else
title(deblank(texname{j}),'interpreter','tex')
end
else
title(deblank(var_list_(j,:)),'interpreter','none')
end
x0=get(gca,'xlim');
hold on, plot(x0, [0 0],'r')
end
if (mod(nplo,9)==0 | j==size(var_list_,1)) & nfigsav<nfig & nfig>0,
nfigsav=nfig;
saveas(gcf,[DirectoryName '\' fname_,'_IRF_comp_',vargin{i},int2str(nfig)])
eval(['print -depsc2 ' DirectoryName '\' fname_,'_IRF_comp_',vargin{i},int2str(nfig)]);
eval(['print -dpdf ' DirectoryName '\' fname_,'_IRF_comp_',vargin{i},int2str(nfig)]);
close(gcf)
end
end
end
nfigsav=0;
nfig=0;
nplo=0;
for j=1:size(var_list_,1),
nplo=nplo+1;
if mod(nplo,9)==1,
figure('name',['Comparison of cumulated orthogonalised shocks'])
nfig=nfig+1;
end
subplot(3,3,nplo-9*(nfig-1))
plot(x,ys(:,1,j),'k',x,ys(:,2,j),':k'), title(deblank(var_list_(j,:)),'interpreter','none')
x0=get(gca,'xlim');
hold on, plot(x0, [0 0],'r')
if (mod(nplo,9)==0 | j==size(var_list_,1)) & nfigsav<nfig & nfig>0,
nfigsav=nfig;
saveas(gcf,[DirectoryName '\' fname_,'_IRF_cum_',int2str(nfig)])
eval(['print -depsc2 ' DirectoryName '\' fname_,'_IRF_cum_',int2str(nfig)]);
eval(['print -dpdf ' DirectoryName '\' fname_,'_IRF_cum_',int2str(nfig)]);
close(gcf)
end
end

82
matlab/dirf_comp_perc.m
View File

@ -1,82 +0,0 @@
function dirf_comp_perc(cname, lname, s1, s2, vargin, texname)
global M_ options_
fname_=M_.fname;
DirectoryName = CheckPath('Output');
if nargin<6,
for j=1:M_.endo_nbr,
texname{j}=deblank(M_.endo_names_tex(j,:));
end
end
for i=1:length(vargin)
nfig=0;
nplo=0;
for j=1:length(cname),
if isfield(s1,[cname{j},'_',vargin{i}]),
nplo=nplo+1;
if mod(nplo,9)==1,
figure('name',['Perc. change comparison of orthogonalised shocks to ',vargin{i}])
nfig=nfig+1;
end
subplot(3,3,nplo-9*(nfig-1))
ytemp_=getfield(s1,[cname{j},'_',vargin{i}]);
y=cumsum(ytemp_);
ytemp_=getfield(s2,[cname{j},'_',vargin{i}]);
y=[y cumsum(ytemp_)];
if ~exist('ys');
ys=zeros([size(y), length(cname)]);
end
ys(:,:,j)=ys(:,:,j)+y;
x=[1:length(y)];
plot(x,y(:,1),'k',x,y(:,2),':k'),
if options_.TeX,
if nargin<6
title(lname{j},'interpreter','none')
else
title(texname{j},'interpreter','tex')
end
else
title(lname{j},'interpreter','none')
end
x0=get(gca,'xlim');
hold on, plot(x0, [0 0],'r')
end
if (mod(nplo,9)==0 | j==length(cname)) & nfig>0,
saveas(gcf,[DirectoryName '\' fname_,'_IRF_comp_perc_',vargin{i},int2str(nfig)])
eval(['print -depsc2 ' DirectoryName '\' fname_,'_IRF_comp_perc_',vargin{i},int2str(nfig)]);
eval(['print -dpdf ' DirectoryName '\' fname_,'_IRF_comp_perc_',vargin{i},int2str(nfig)]);
close(gcf)
end
end
end
nfigsav=0;
nfig=0;
nplo=0;
for j=1:length(cname),
nplo=nplo+1;
if mod(nplo,9)==1,
figure('name',['Perc. comparison of cumulated orthogonalised shocks'])
nfig=nfig+1;
end
subplot(3,3,nplo-9*(nfig-1))
plot(x,ys(:,1,j),'k',x,ys(:,2,j),':k'),
if options_.TeX,
title(texname{j},'interpreter','tex')
else
title(lname{j},'interpreter','none')
end
x0=get(gca,'xlim');
hold on, plot(x0, [0 0],'r')
if (mod(nplo,9)==0 | j==length(cname)) & nfigsav<nfig & nfig>0,
nfigsav=nfig;
saveas(gcf,[DirectoryName '\' fname_,'_IRF_cum_perc_',int2str(nfig)])
eval(['print -depsc2 ' DirectoryName '\' fname_,'_IRF_cum_perc_',int2str(nfig)]);
eval(['print -dpdf ' DirectoryName '\' fname_,'_IRF_cum_perc_',int2str(nfig)]);
close(gcf)
end
end

57
matlab/dirf_comp_perc_post.m
View File

@ -1,57 +0,0 @@
function dirf_comp_perc_post(lname, s1, s2, varargin)
global M_ options_
fname_=M_.fname;
DirectoryName = CheckPath('Output');
in=size(lname);
if min(in)==2,
texname=lname(2,:);
lname=lname(1,:);
else
for j=1:length(lname),
iendo=strmatch(lname{j},M_.endo_names,'exact');
texname{j}=deblank(M_.endo_names_tex(iendo,:));
end
end
for i=1:length(varargin)
nfig=0;
nplo=0;
for j=1:length(lname),
y=getfield(s1.PosteriorIRF.Mean,[lname{j},'_',varargin{i}]);
y=[y getfield(s2.PosteriorIRF.Mean,[lname{j},'_',varargin{i}])];
if max(max(abs(y)))>1e-10 & max(max(y)-min(y))>1e-10;
nplo=nplo+1;
if mod(nplo,9)==1,
figure('name',['Perc. change comparison of orthogonalised shocks to ',varargin{i}])
nfig=nfig+1;
end
subplot(3,3,nplo)
yl=getfield(s1.PosteriorIRF.HPDinf,[lname{j},'_',varargin{i}]);
yl=[yl getfield(s2.PosteriorIRF.HPDinf,[lname{j},'_',varargin{i}])];
yu=getfield(s1.PosteriorIRF.HPDsup,[lname{j},'_',varargin{i}]);
yu=[yu getfield(s2.PosteriorIRF.HPDsup,[lname{j},'_',varargin{i}])];
x=[1:length(y)];
% patch([x, x(end:-1:1)], [yu(: ,1)', yl(end:-1:1,1)'],[0 0.75 0],'facealpha',0.5);
% hold on,
% patch([x, x(end:-1:1)], [yu(: ,2)', yl(end:-1:1,2)'],[0.75 0 0],'facealpha',0.5);
plot(x,y(:,1),'k',x,y(:,2),':k'),
hold on,
if options_.TeX,
title(texname{j},'interpreter','tex')
else
title(lname{j},'interpreter','none')
end
plot(x([1,end]), [0 0],'r')
hold off,
end
if (mod(nplo,9)==0 | j==length(lname)) & nfig>0,
saveas(gcf,[DirectoryName,'/',fname_,'_Bayesian_IRF_comp_perc_',varargin{i},'_',int2str(nfig)])
eval(['print -depsc2 ' DirectoryName,'/',fname_,'_Bayesian_IRF_comp_perc_',varargin{i},'_',int2str(nfig)]);
eval(['print -dpdf ' DirectoryName,'/',fname_,'_Bayesian_IRF_comp_perc_',varargin{i},'_',int2str(nfig)]);
close(gcf)
nplo=0;
end
end
end

49
matlab/dirf_comp_post.m
View File

@ -1,49 +0,0 @@
function dirf_comp_post(var_list_, s1, s2, vargin, texname)
global M_ options_
fname_ = M_.fname;
DirectoryName = CheckPath('Output');
for i=1:length(vargin)
nfigsav=0;
nfig=0;
nplo=0;
for j=1:length(var_list_),
y=getfield(s1.PosteriorIRF.Mean,[var_list_{j},'_',vargin{i}]);
y=[y getfield(s2.PosteriorIRF.Mean,[var_list_{j},'_',vargin{i}])];
if max(max(abs(y)))>1.e-10,
nplo=nplo+1;
if mod(nplo,9)==1,
figure('name',['Comparison of orthogonalised shocks to ',vargin{i}])
nfig=nfig+1;
end
subplot(3,3,nplo)
yl=getfield(s1.PosteriorIRF.HPDinf,[var_list_{j},'_',vargin{i}]);
yu=getfield(s1.PosteriorIRF.HPDsup,[var_list_{j},'_',vargin{i}]);
yl=[yl getfield(s2.PosteriorIRF.HPDinf,[var_list_{j},'_',vargin{i}])];
yu=[yu getfield(s2.PosteriorIRF.HPDsup,[var_list_{j},'_',vargin{i}])];
x=[1:size(y,1)];
patch([x, x(end:-1:1)], [yu(: ,1)', yl(end:-1:1,1)'],[0.75 0.75 0.75]); %,'facealpha',0.5);
hold on,
%patch([x, x(end:-1:1)], [yu(: ,2)', yl(end:-1:1,2)'],[0.75 0 0],'facealpha',0.5);
plot(x,y(:,1),'k',x,y(:,2),'--k','LineWidth',1),
if options_.TeX & nargin==5
title(texname{j},'interpreter','tex')
else
title(var_list_{j},'interpreter','none')
end
plot(x([1, end]), [0 0],'r')
hold off,
end
if (mod(nplo,9)==0 | j==length(var_list_)) & nfigsav<nfig & nfig>0,
nfigsav=nfig;
saveas(gcf,[DirectoryName,'\',fname_,'_Bayesian_IRF_comp_',vargin{i},'_',int2str(nfig)])
eval(['print -depsc2 ' DirectoryName,'\',fname_,'_Bayesian_IRF_comp_',vargin{i},'_',int2str(nfig)]);
eval(['print -dpdf ' DirectoryName,'\',fname_,'_Bayesian_IRF_comp_',vargin{i},'_',int2str(nfig)]);
close(gcf)
nplo=0;
end
end
end

137
matlab/dlr01.m
View File

@ -1,137 +0,0 @@
function [fit,SE,par,V,COMP,o1,y0]= dlr0(y,z,IRW,nvr,alpha,P0,x0,smooth,ALG,p)
% Version for use with SDP01 for fast single parameter estimation
% DLR Dynamic Linear Regression analysis
%
% [fit,fitse,par,parse,comp,e,y0]=dlr(y,z,TVP,nvr,alpha,P0,x0,sm,ALG)
% 1 2 3 4 5 6 7 8 9
%
% y: Time series (*)
% z: Regressors (*)
% TVP: Model type for each TVP (0-RW/AR(1), 1-IRW/SRW) (0)
% nvr: NVR hyper-parameters (0)
% alpha: alpha hyper-parameters for SRW model (1)
% P0: Initial P matrix (1e5)
% x0: Initial state vector (0)
% sm: Smoothing on (1-default) or off (0-saves memory)
% ALG: Smoothing algorithm: P (0) or Q (1-default)
%
% fit: Model fit
% fitse: Standard error of the fit
% par: Parameter estimates
% parse: Standard errors of parameters
% comp: Linear components
% e: Normalised innovations; use e=e(~isnan(e)) to remove NaNs
% y0: Interpolated data
%
% Example: dlr(y, [ones(size(u)) u], [0 1], 0.001, [0.95 1])
% regression type model y = c1 + c2(t)*u, with an AR(1) model for
% c1 (alpha=0.95) and an IRW model for c2 (NVR=0.001 in both cases)
%
% See also DLROPT, FCAST, STAND, DHR, DHROPT, SDP
% Copyright (c) 2004 by CRES, Lancaster University, United Kingdom
% Authors : Peter Young, Wlodek Tych, Diego Pedregal, James Taylor
% Revision history :
% 04/03/2003, JT, cosmetic
% 22/08/2002, JT, interpolated data returned
% 18/01/1999, JT, revised output arguments
% 02/11/1998, JT, fit output argument
% 23/09/1998, JT, separate y and regressors input arguments, Q/P algorithm
% 22/09/1998, JT, SRW option, initial P0 not overwritten
% 01/10/1997, DP, univariate toolbox
% Auxiliary functions:
% dlr_opt:
% dlrfun:
y=y(:); % ensure y is a column, 28/6/99, JT
Z=[y, z]; %23/9/98, JT
[n, m]= size(Z);
% Checking inputs
if nargin<10, p=[]; end
if nargin<9, ALG= 1; end %23/9/98, JT
if isempty(ALG), ALG= 1; end
if nargin<8, smooth= 1; end
if isempty(smooth), smooth=1; end
if nargin<7, x0= []; end
if nargin<6, P0= []; end
if ~isempty(P0), lP0= length(P0);
P0= P0(:)'; P0= [P0 ones(1, m-1-lP0)*P0(lP0)];
end
if nargin<5, alpha= 1; end %JT, 22/9/98, SRW option
if isempty(alpha), alpha= 1; end
if any(size(P0))==1
P0= diag(P0);
end
% Checking missing values
I= zeros(n, 1); maxp= m-1;
Miss= I;
% SS F, Q and P0 matrices
if IRW
F=[1 1 ; 0 1]; Q=[ 0 0 ; 0 nvr];
else
F=1; Q=nvr;
end
[row, col]= size(F);
if isempty(P0) %JT, 22/09/98, don't overwrite input argument
P0= eye(row)*1e5;
end
% Filtering and Smoothing
[Xkk, Pkk, er, ykk1, Xkk1, yp, ep, ykk, XkN, ys, PkN, ers]=ksirw2c(Z, IRW, Miss, I, F, P0, Q, 0, x0, [], [], ALG);
Xkk= XkN; Pkk= PkN; er= ers;
% Building outputs (linear components, parameters, ...)
i= cumsum([1 IRW(1:maxp-1)+1]);
% TVP parameters
par= Xkk(i, :)';
% Linear components
COMP= par.*Z(:, 2:m);
% Normalised innovations and SE of predictions
o1= [ones(row,1)*nan; (Z(row+1:n, 1)-ykk1(row+1:n))./sqrt(er(row+1:n))]; % Normalised innovations
o11= o1(~isnan(o1));
sigma= o11'*o11/(length(o11)); % Innovations variance
SE= sqrt(sigma.*er); % SE of innovations
% Standard error of parameters
Pkk= Pkk';
W= Pkk*sigma;
V= zeros(n, m-1); [rP, cP]= size(Pkk);
for j= 1:length(i)
V(:, j)= sqrt(W(i(j):cP:rP, i(j)));
end
% Fit, 2/10/98, JT
if size(COMP, 2)>1
fit=sum(COMP')';
else
fit=COMP;
end
% interpolated data, 22/08/02, JT
ii=find(isnan(y));
y0=y;
y0(ii)=fit(ii);
% end of m-file

BIN
matlab/dlr01.p
View File

Binary file not shown.

79
matlab/dlrfun01.m
View File

@ -1,79 +0,0 @@
function [out,er,o1]= dlrfun(b0, X, nvr, RWtype, Miss, F, P0, steps, ALG, p, stop_info)
%DLRFUN Likelihood function for DLR hyper-parameter estimation.
% out= dlrfun(b0, X, nvr, RWtype, Miss, F, P0, steps, ALG, p, stop_info)
%
% out: Value of the likelihood function.
%
% b0: Alphas and/or NVR (transformed) parameters for evaluation of the likelihood.
% X: Time series (output, inputs).
% nvr: Constrained estimation (number of inputs vector). Options:
% -2: Free estimation of nvr parameters (default).
% -1: Constrained estimation (all the nvrs at locations
% where nvr is -1 will be the same).
% any value > 0: Nvr constrained to such value (not estimated).
% RWtype: RW type of TVP parameter (0-RW, 1-IRW).
% Miss: Location of missing values
% F: SS transition matrix. When SRW models required, the diagonal of F may contain:
% -2: Free estimation of alpha parameters.
% -1: Constrained estimation (all the alphas at locations
% where F(i, i) is -1 will be the same).
% any value > 0: Alpha constrained to such value (default: 1).
% P0: Initial P matrix.
% steps: Number of steps ahead forecasts for objective (0: ML)
% ALG: Optimisation algorithm: 0=fmins, 1=fminu, 2=leastsq (not ML) (0).
% p: extra parameter for DAR modelling
% stop_info: graphics handle for optimisation window
% Copyright (c) 2004 by CRES, Lancaster University, United Kingdom
% Authors : Peter Young, Wlodek Tych, Diego Pedregal, James Taylor
[n, cX]= size(X);
NVR= 10.^(b0);
% Constructing nvr vector
x= NVR;
i1= min(find(nvr== -1)); i2= sum(nvr(1:i1)== -2);
iii= find(nvr<0); jjj= find(nvr>-1);
if ~isempty(i1), rest= i2+1; else, rest= 0; end
j= 1;
for i= 1:cX-1
if nvr(i)== -2
if j== rest, j= j+1; end
nvr(i)= x(j); j= j+1;
end
if nvr(i)== -1, nvr(i)= x(rest); end
end
if ~isempty(iii), x(iii)= 10.^(nvr(iii)); end
if ~isempty(jjj), x(jjj)= nvr(jjj); end
x= x(:);
if RWtype;
Q=[0 0;0 NVR];
else
Q=NVR;
end
[Xkk, Pkk, er, ykk1, Xkk1, yp, ep, ykk]= ...
ksirw2c(X, RWtype, Miss, zeros(size(Miss)), F, P0, Q, 0, [], steps);
ind= ((length(F)+2):length(er))';
nd=length(ind);
out= (X(ind, 1)-ykk1(ind))./sqrt(er(ind));
out= sum(log(er(ind)))+nd*log((out'*out)/nd);
if nargout>1, o1= X(ind, 1)-ykk1(ind); er=er(ind); end
%JT, 9/4/99
if stop_info
if steps==0 % ML Log-Likelihood
set(stop_info, 'String', num2str(-0.5*out-size(X, 1)/2*(log(2*pi)+1)))
elseif ALG==2
set(stop_info, 'String', num2str(out'*out)) % display scalar output
else
set(stop_info, 'String', num2str(out))
end
drawnow
end
% end of m-file

BIN
matlab/dlrfun01.p
View File

Binary file not shown.

141
matlab/dlropt01.m
View File

@ -1,141 +0,0 @@
function [NVR,ALPHA,opts,separ]= dlropt0(y,z,IRW,method,nvr,alpha,nvr0,alpha0,opts,ALG,output,P0,p);
% DLROPT Hyper-parameter estimation for DLR
%
% [nvr,alpha,opts,parse]=dlropt(y,z,TVP,meth,nvrc,alphac,nvr0,alpha0,opts,ALG,tab,P0)
% 1 2 3 4 5 6 7 8 9 10 11 12
%
% y: Time series (*)
% z: Regressors (*)
% TVP: Model type for each TVP (0-RW/AR(1), 1-IRW/SRW) (0)
% meth: Estimation method ('ml')
% 'ml': Maximum Likelihood
% 'f#': Sum of squares of the #-step-ahead forecasting errors
% nvrc: Constraints for each NVR (-2)
% -2: Free estimation
% -1: Constrained estimation (all parameters with nvrc=-1 are equal)
% >=0: NVR constrained to this value (it is not estimated)
% alphac: Constraints for each alpha (-2, -1, or >=0 as for nvrc) (1)
% nvr0: Initial NVR hyper-parameters (0.0001)
% alpha0: Initial alpha hyper-parameters (1)
% opts: Optimisation options. Type 'help foptions' for details
% ALG: Optimisation algorithm: 0=fmins, 1=fminu, 2=leastsq (not ML) (0)
% tab: Display: 0=none, 1=tabulate results, 2=update window (2)
% P0: Initial P matrix (1e5)
%
% nvr: Estimated NVR hyper-parameters
% alpha: Estimated alpha hyper-parameters
% opts: Returned optimisation options. Type 'help foptions' for details
% parse: Standard Error of hyper-parameters (omit to reduce computation time)
%
% Example: dlropt(y, [ones(size(u)) u], 0, [], [0 -2])
% regression type model y = c1 + c2(t)*u, with an RW model for
% both parameters and c1 assumed constant (NVR fixed at zero)
%
% See also DLR, FCAST, STAND, FOPTIONS, FMINS
% Copyright (c) 2004 by CRES, Lancaster University, United Kingdom
% Authors : Peter Young, Wlodek Tych, Diego Pedregal, James Taylor
y=y(:); % ensure y is a column, 28/6/99, JT
Z=[y, z]; %23/9/98, JT
[n, cZ]= size(Z);
p=[];
steps=0;
P0= diag(P0);
zlabel='Z1';
nnvr=1;
Miss= isnan(Z(:,1));
F= zeros(sum(IRW+1)); j=1;
for i= 1:cZ-1
if IRW(i), f= [alpha(i) 1; 0 1]; else, f= alpha(i); end
F(j:j+IRW(i), j:j+IRW(i))= f; j= j+IRW(i)+1;
end
% Initial conditions ... strange initial condition for no alpha
% optimisation .... ******************* WT 11/04
%b0= log10([alpha0./(1-alpha0) nvr0]);
b0=log10(nvr0);
t1= clock;
% Optimisation update windows, JT, revised 3/8/99
% ALG is 0
% Matlab is newer than 5
% output is 2
% steps is 0
if output==2
update_panel=figure('Units', 'Normalized','Position', [0.425 0.5 0.15 0.15], ...
'Number', 'Off', 'Name', 'Working...', 'MenuBar', 'None');
stop_alg=uicontrol(update_panel, 'Style', 'Text', ...
'Units', 'Normalized', 'Position', [0.1 12/15 0.8 2/15]);
stop_meth=uicontrol(update_panel, 'Style', 'Text', ...
'Units', 'Normalized', 'Position', [0.1 9/15 0.8 2/15]);
stop_info=uicontrol(update_panel, 'Style', 'Text', 'Units', ...
'Normalized', 'Position', [0.1 6/15 0.8 2/15], 'Tag', 'info');
stop_but=uicontrol(update_panel, 'Style', 'Push', 'String', 'S T O P', ...
'Call', 'set(gcbf,''UserData'', 1);', 'Units', 'Normalized', ...
'Position', [0.1 1/15 0.8 4/15]);
set(update_panel, 'UserData', 0); % before stop button is pressed
set(stop_alg, 'String', 'ALG = fmins')
set(stop_meth, 'String', 'Log-likelihood')
end
if exist('stop_info')==0, stop_info=[]; end
[lb, opts]= fstop('dlrfun01', b0, opts, [], ...
Z, nvr, IRW, Miss, F, P0, steps, ALG, p, stop_info); % JT, 30/01/04, fstop
if exist('update_panel'), close(update_panel), end
separ= nan*ones(size(lb));
separ=seml('dlrfun01', lb, Z, nvr, IRW, Miss, F, P0, steps, ALG, p, 0);
tiempo= etime(clock, t1);
NVR=lb;
SEPARn=separ;
SEPARa=0;
ALPHA= alpha';
NVR=10.^NVR;
% Display results
if output
disp(' ')
if tiempo>60
t1= [int2str(fix(tiempo/60)) ' minutes, ' ...
num2str(fix(rem(tiempo, 60))) ' seconds.'];
else, t1= [num2str(tiempo) ' seconds.'];
end
hora= fix(clock); if hora(5)<10, pega= '0'; else, pega= []; end
disp('TIME DOMAIN ESTIMATION')
if steps==0, disp('METHOD: MAXIMUM LIKELIHOOD'), % DP, 10/02/99
else, disp(['METHOD: SUM OF SQUARES ' int2str(steps) '-STEPS-AHEAD FORECAST ERRORS'])
end
if ALG==0, disp('OPTIMISER: FMINS'); % DP, 10/02/99
elseif ALG==1, disp('OPTIMISER: FMINU');
elseif ALG==2, disp('OPTIMISER: LEASTSQ');
end
disp(['Date: ' date ' / Time: ' num2str(hora(4)) ':' pega num2str(hora(5))])
disp(t1)
disp('RW NVR -S.E. +S.E.')
disp([num2str(IRW,1) ' ' num2str(NVR,4) ' ' ,...
num2str(10.^(lb-SEPARn),3) ' ' num2str(10.^(lb+SEPARn),3)])
if steps==0, disp(['Log-likelihood: ' num2str(-0.5*opts(8)-length(y)/2*(log(2*pi)+1))])
else, disp(['Sum of Squares of ' int2str(steps) '-steps-ahead forecast errors: ' num2str(opts(8))])
end
disp(' ')
end
% end of m-file

BIN
matlab/dlropt01.p
View File

Binary file not shown.

202
matlab/dyn_bvar.m
View File

@ -1,202 +0,0 @@
function dyn_bvar(maxnlags,tau,d,lambda,mu,omega,train,flat,breaks)
% stephane.adjemian@cepremap.cnrs.fr [21 april, 2004]
global options_
eval(options_.datafile);
dataset = [ ];
for i=1:size(options_.varobs,1)
dataset = [dataset eval(deblank(options_.varobs(i,:)))];
end
default_nlags = 8;
default_const = 1;
default_tau = 3;
default_d = 0.5;
default_lambda = 5;
default_mu = 2;
default_omega = 1;
if nargin == 0
maxnlags = default_nlags;
tau = default_tau;
d = default_d;
lambda = default_lambda;
mu = default_mu;
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 1
tau = default_tau;
d = default_d;
lambda = default_lambda;
mu = default_mu;
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 2
if isempty(maxnlags)
maxnlags = default_nlags;
end
d = default_d;
lambda = default_lambda;
mu = default_mu;
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 3
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(tau)
tau = default_tau;
end
lambda = default_lambda;
mu = default_mu;
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 4
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
mu = default_mu;
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 5
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
if isempty(lambda)
lambda = default_lambda;
end
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 6
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
if isempty(lambda)
lambda = default_lambda;
end
if isempty(mu)
mu = default_mu;
end
train = [];
flat = 0;
breaks = [];
elseif nargin == 7
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
if isempty(lambda)
lambda = default_lambda;
end
if isempty(mu)
mu = default_mu;
end
if isempty(omega)
omega = default_omega;
end
flat = 0;
breaks = [];
elseif nargin == 8
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(const)
const = default_const;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
if isempty(lambda)
lambda = default_lambda;
end
if isempty(mu)
mu = default_mu;
end
if isempty(omega)
omega = default_omega;
end
breaks = [];
train = [];
elseif nargin == 9
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(const)
const = default_const;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
if isempty(lambda)
lambda = default_lambda;
end
if isempty(mu)
mu = default_mu;
end
if isempty(omega)
omega = default_omega;
end
train = [];
flat = 0;
elseif nargin > 9
disp('dyn_bvar :: too many arguments.')
end
mnprior.tight = tau;
mnprior.decay = d;
vprior.sig = std(dataset(options_.first_obs+options_.presample-maxnlags:options_.first_obs+options_.presample-1,:))';
vprior.w = omega;
for lag = 1:maxnlags
ydata = dataset(options_.first_obs+options_.presample-lag:options_.first_obs+options_.nobs-1,:);
w=mgnldnsty(ydata,lag,ones(options_.nobs+lag-options_.presample,1),breaks,lambda,mu,mnprior,vprior,train,flat);
disp(' ')
fprintf('The marginal log density of the BVAR(%g) model is equal to %10.4f \n',lag,w);
disp(' ')
end

207
matlab/dyn_bvar2.m
View File

@ -1,207 +0,0 @@
function dyn_bvar(maxnlags,tau,d,lambda,mu,omega,train,flat,breaks)
% stephane.adjemian@cepremap.cnrs.fr [22 april, 2004]
global options_
eval(options_.datafile);
dataset = [ ];
for i=1:size(options_.varobs,1)
dataset = [dataset eval(deblank(options_.varobs(i,:)))];
end
default_nlags = 8;
default_const = 1;
default_tau = 3;
default_d = 0.5;
default_lambda = 5;
default_mu = 2;
default_omega = 1;
if nargin == 0
maxnlags = default_nlags;
tau = default_tau;
d = default_d;
lambda = default_lambda;
mu = default_mu;
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 1
tau = default_tau;
d = default_d;
lambda = default_lambda;
mu = default_mu;
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 2
if isempty(maxnlags)
maxnlags = default_nlags;
end
d = default_d;
lambda = default_lambda;
mu = default_mu;
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 3
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(tau)
tau = default_tau;
end
lambda = default_lambda;
mu = default_mu;
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 4
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
mu = default_mu;
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 5
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
if isempty(lambda)
lambda = default_lambda;
end
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 6
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
if isempty(lambda)
lambda = default_lambda;
end
if isempty(mu)
mu = default_mu;
end
train = [];
flat = 0;
breaks = [];
elseif nargin == 7
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
if isempty(lambda)
lambda = default_lambda;
end
if isempty(mu)
mu = default_mu;
end
if isempty(omega)
omega = default_omega;
end
flat = 0;
breaks = [];
elseif nargin == 8
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(const)
const = default_const;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
if isempty(lambda)
lambda = default_lambda;
end
if isempty(mu)
mu = default_mu;
end
if isempty(omega)
omega = default_omega;
end
breaks = [];
train = [];
elseif nargin == 9
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(const)
const = default_const;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
if isempty(lambda)
lambda = default_lambda;
end
if isempty(mu)
mu = default_mu;
end
if isempty(omega)
omega = default_omega;
end
train = [];
flat = 0;
elseif nargin > 9
disp('dyn_bvar :: too many arguments.')
end
mnprior.tight = tau;
mnprior.decay = d;
yy = dataset(1:options_.first_obs+options_.presample-1,:);
bb = inv(yy(1:size(yy,1)-1,:)'*yy(1:size(yy,1)-1,:))*yy(1:size(yy,1)-1,:)'*yy(2:size(yy,1),:);
rr = yy(2:size(yy,1),:)-yy(1:size(yy,1)-1,:)*bb;
vprior.sig = sqrt(diag(rr'*rr/(size(yy,1)-1)))';
vprior.w = omega;
% mgnldnsty(ydata,lags,xdata,breaks,lambda,mu,mnprior,vprior,train,flat)
for lag = 1:maxnlags
%ydata = dataset(options_.first_obs+options_.presample-lag:options_.first_obs+options_.presample+options_.nobs-1,:);
ydata = dataset(options_.first_obs+options_.presample-lag:options_.first_obs+options_.nobs-1,:);
w=mgnldnsty(ydata,lag,ones(options_.nobs+lag-options_.presample,1),breaks,lambda,mu,mnprior,vprior,train,flat);
disp(' ')
fprintf('The marginal log density of the BVAR(%g) model is equal to %10.4f \n',lag,w);
disp(' ')
end

207
matlab/dyn_bvar_mr.m
View File

@ -1,207 +0,0 @@
function dyn_bvar(maxnlags,tau,d,lambda,mu,omega,train,flat,breaks)
% stephane.adjemian@cepremap.cnrs.fr [21 april, 2004]
global options_
eval(options_.datafile);
dataset = [ ];
for i=1:size(options_.varobs,1)
dataset = [dataset eval(deblank(options_.varobs(i,:)))];
end
default_nlags = 8;
default_const = 1;
default_tau = 3;
default_d = 0.5;
default_lambda = 5;
default_mu = 2;
default_omega = 1;
if nargin == 0
maxnlags = default_nlags;
tau = default_tau;
d = default_d;
lambda = default_lambda;
mu = default_mu;
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 1
tau = default_tau;
d = default_d;
lambda = default_lambda;
mu = default_mu;
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 2
if isempty(maxnlags)
maxnlags = default_nlags;
end
d = default_d;
lambda = default_lambda;
mu = default_mu;
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 3
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(tau)
tau = default_tau;
end
lambda = default_lambda;
mu = default_mu;
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 4
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
mu = default_mu;
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 5
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
if isempty(lambda)
lambda = default_lambda;
end
omega = default_omega;
breaks = [];
train = [];
flat = 0;
elseif nargin == 6
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
if isempty(lambda)
lambda = default_lambda;
end
if isempty(mu)
mu = default_mu;
end
train = [];
flat = 0;
breaks = [];
elseif nargin == 7
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
if isempty(lambda)
lambda = default_lambda;
end
if isempty(mu)
mu = default_mu;
end
if isempty(omega)
omega = default_omega;
end
flat = 0;
breaks = [];
elseif nargin == 8
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(const)
const = default_const;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
if isempty(lambda)
lambda = default_lambda;
end
if isempty(mu)
mu = default_mu;
end
if isempty(omega)
omega = default_omega;
end
breaks = [];
train = [];
elseif nargin == 9
if isempty(maxnlags)
maxnlags = default_nlags;
end
if isempty(const)
const = default_const;
end
if isempty(tau)
tau = default_tau;
end
if isempty(d)
d = default_d;
end
if isempty(lambda)
lambda = default_lambda;
end
if isempty(mu)
mu = default_mu;
end
if isempty(omega)
omega = default_omega;
end
train = [];
flat = 0;
elseif nargin > 9
disp('dyn_bvar :: too many arguments.')
end
%load cul;
%dataset=xx;
% tried with VECM!!!!
mnprior.tight = tau;
mnprior.decay = d;
vprior.sig = std(dataset(options_.first_obs+options_.presample-maxnlags:options_.first_obs+options_.presample-1,:))';
vprior.w = omega;
xdata=[];
for lag = 1:maxnlags
ydata = dataset(options_.first_obs+options_.presample-lag:options_.first_obs+options_.nobs-1,:);
xdata = xvec(options_.first_obs+options_.presample-lag:options_.first_obs+options_.presample+options_.nobs-1,:);
w=mgnldnsty(ydata,lag,[ones(options_.nobs+lag,1) xdata],breaks,lambda,mu,mnprior,vprior,train,flat);
disp(' ')
fprintf('The marginal log density of the BVAR(%g) model is equal to %10.4f \n',lag,w);
disp(' ')
end

309
matlab/dynare_MC.m
View File

@ -1,309 +0,0 @@
function dynare_MC(var_list_,OutDir)
global M_ options_ oo_ estim_params_
global bayestopt_
% temporary fix until M_.H is initialized by the parser
M_.H = [];
options_.varlist = var_list_;
options_.lgyidx2varobs = zeros(size(M_.endo_names,1),1);
for i = 1:size(M_.endo_names,1)
tmp = strmatch(deblank(M_.endo_names(i,:)),options_.varobs,'exact');
if ~isempty(tmp)
options_.lgyidx2varobs(i,1) = tmp;
end
end
options_ = set_default_option(options_,'first_obs',1);
options_ = set_default_option(options_,'prefilter',0);
options_ = set_default_option(options_,'presample',0);
options_ = set_default_option(options_,'lik_algo',1);
options_ = set_default_option(options_,'lik_init',1);
options_ = set_default_option(options_,'nograph',0);
options_ = set_default_option(options_,'mh_conf_sig',0.90);
options_ = set_default_option(options_,'mh_replic',20000);
options_ = set_default_option(options_,'mh_drop',0.5);
options_ = set_default_option(options_,'mh_jscale',0.2);
options_ = set_default_option(options_,'mh_init_scale',2*options_.mh_jscale);
options_ = set_default_option(options_,'mode_file','');
options_ = set_default_option(options_,'mode_compute',4);
options_ = set_default_option(options_,'mode_check',0);
options_ = set_default_option(options_,'prior_trunc',1e-10);
options_ = set_default_option(options_,'mh_mode',1);
options_ = set_default_option(options_,'mh_nblck',2);
options_ = set_default_option(options_,'load_mh_file',0);
options_ = set_default_option(options_,'nodiagnostic',0);
options_ = set_default_option(options_,'loglinear',0);
options_ = set_default_option(options_,'unit_root_vars',[]);
options_ = set_default_option(options_,'XTick',[]);
options_ = set_default_option(options_,'XTickLabel',[]);
options_ = set_default_option(options_,'bayesian_irf',0);
options_ = set_default_option(options_,'bayesian_th_moments',0);
options_ = set_default_option(options_,'TeX',0);
options_ = set_default_option(options_,'irf',40);
options_ = set_default_option(options_,'relative_irf',0);
options_ = set_default_option(options_,'order',1);
options_ = set_default_option(options_,'ar',5);
options_ = set_default_option(options_,'dr_algo',0);
options_ = set_default_option(options_,'linear',0);
options_ = set_default_option(options_,'drop',0);
options_ = set_default_option(options_,'replic',1);
options_ = set_default_option(options_,'hp_filter',0);
options_ = set_default_option(options_,'forecast',0);
options_ = set_default_option(options_,'smoother',0);
options_ = set_default_option(options_,'moments_varendo',0);
options_ = set_default_option(options_,'filtered_vars',0);
options_ = set_default_option(options_,'kalman_algo',1);
options_ = set_default_option(options_,'kalman_tol',10^(-12));
options_ = set_default_option(options_,'posterior_mode_estimation',1);
options_ = set_default_option(options_,'MaxNumberOfBytes',1e6);
options_ = set_default_option(options_,'xls_sheet','');
options_ = set_default_option(options_,'xls_range','');
options_ = set_default_option(options_,'filter_step_ahead',0);
options_ = set_default_option(options_,'diffuse_d',[]);
options_ = set_default_option(options_,'Opt6Iter',3);
options_ = set_default_option(options_,'Opt6Numb',100000);
options_ = set_default_option(options_,'steadystate_flag',0);
options_ = set_default_option(options_,'logdata',0);
options_ = set_default_option(options_,'use_mh_covariance_matrix',0);
if exist([M_.fname '_steadystate.m'])
options_.steadystate_flag = 1;
end
if options_.filtered_vars ~= 0 & options_.filter_step_ahead == 0
options_.filter_step_ahead = 1;
end
if options_.bayesian_irf & options_.irf == 0
options_.irf = 40;
end
if options_.filter_step_ahead ~= 0
options_.nk = max(options_.filter_step_ahead);
else
options_.nk = 0;
end
%% Add something to the parser ++>
M_.dname = M_.fname; % The user should be able to choose another name
% for the directory...
pnames = [' ';'beta ';'gamm ';'norm ';'invg ';'unif ';'invg2'];
n_varobs = size(options_.varobs,1);
[xparam1,estim_params_,bayestopt_,lb,ub] = set_prior(estim_params_);
bounds = prior_bounds(bayestopt_);
bounds(:,1)=max(bounds(:,1),lb);
bounds(:,2)=min(bounds(:,2),ub);
if any(xparam1 < bounds(:,1)) | any(xparam1 > bounds(:,2))
find(xparam1 < bounds(:,1))
find(xparam1 > bounds(:,2))
error('Initial parameter values are outside parameter bounds')
end
lb = bounds(:,1);
ub = bounds(:,2);
bayestopt_.lb = lb;
bayestopt_.ub = ub;
if ~isfield(options_,'trend_coeffs')
bayestopt_.with_trend = 0;
else
bayestopt_.with_trend = 1;
bayestopt_.trend_coeff = {};
trend_coeffs = options_.trend_coeffs;
nt = length(trend_coeffs);
for i=1:n_varobs
if i > length(trend_coeffs)
bayestopt_.trend_coeff{i} = '0';
else
bayestopt_.trend_coeff{i} = trend_coeffs{i};
end
end
end
bayestopt_.penalty = 1e8; % penalty
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
nx = nvx+nvn+ncx+ncn+np;
dr = set_state_space([]);
%% Initialization with unit-root variables
if ~isempty(options_.unit_root_vars)
n_ur = size(options_.unit_root_vars,1);
i_ur = zeros(n_ur,1);
for i=1:n_ur
i1 = strmatch(deblank(options_.unit_root_vars(i,:)),M_.endo_names(dr.order_var,:),'exact');
if isempty(i1)
error('Undeclared variable in unit_root_vars statement')
end
i_ur(i) = i1;
end
if M_.maximum_lag > 1
l1 = flipud([cumsum(M_.lead_lag_incidence(1:M_.maximum_lag-1,dr.order_var),1);ones(1,M_.endo_nbr)]);
n1 = nnz(l1);
bayestopt_.Pinf = zeros(n1,n1);
l2 = find(l1');
l3 = zeros(M_.endo_nbr,M_.maximum_lag);
l3(i_ur,:) = l1(:,i_ur)';
l3 = l3(:);
i_ur1 = find(l3(l2));
i_stable = ones(M_.endo_nbr,1);
i_stable(i_ur) = zeros(n_ur,1);
i_stable = find(i_stable);
bayestopt_.Pinf(i_ur1,i_ur1) = diag(ones(1,length(i_ur1)));
bayestopt_.i_var_stable = i_stable;
l3 = zeros(M_.endo_nbr,M_.maximum_lag);
l3(i_stable,:) = l1(:,i_stable)';
l3 = l3(:);
bayestopt_.i_T_var_stable = find(l3(l2));
else
n1 = M_.endo_nbr;
bayestopt_.Pinf = zeros(n1,n1);
bayestopt_.Pinf(i_ur,i_ur) = diag(ones(1,length(i_ur)));
l1 = ones(M_.endo_nbr,1);
l1(i_ur,:) = zeros(length(i_ur),1);
bayestopt_.i_T_var_stable = find(l1);
end
options_.lik_init = 3;
end % if ~isempty(options_.unit_root_vars)
if isempty(options_.datafile)
error('ESTIMATION: datafile option is missing')
end
if isempty(options_.varobs)
error('ESTIMATION: VAROBS is missing')
end
%% If jscale isn't specified for an estimated parameter, use
%% global option options_.jscale, set to 0.2, by default
k = find(isnan(bayestopt_.jscale));
bayestopt_.jscale(k) = options_.mh_jscale;
%% Read and demean data
rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
k = [];
k1 = [];
for i=1:n_varobs
k = [k strmatch(deblank(options_.varobs(i,:)),M_.endo_names(dr.order_var,:),'exact')];
k1 = [k1 strmatch(deblank(options_.varobs(i,:)),M_.endo_names, 'exact')];
end
% union of observed and state variables
k2 = union(k',[dr.nstatic+1:dr.nstatic+dr.npred]');
% including variables in t-2 and earlier, if any
k2 = [k2;[M_.endo_nbr+(1:dr.nspred-dr.npred)]'];
% set restrict_state to postion of observed + state variables
% in expanded state vector
bayestopt_.restrict_state = k2;
% set mf1 to positions of observed variables in restricted state vector
% for likelihood computation
[junk,bayestopt_.mf1] = ismember(k,k2);
% set mf2 to positions of observed variables in expanded state vector
% for filtering and smoothing
bayestopt_.mf2 = k;
bayestopt_.mfys = k1;
options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1);
gend = options_.nobs;
rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
if options_.loglinear == 1 & ~options_.logdata
rawdata = log(rawdata);
end
if options_.prefilter == 1
bayestopt_.mean_varobs = mean(rawdata,1);
data = transpose(rawdata-ones(gend,1)*bayestopt_.mean_varobs);
else
data = transpose(rawdata);
end
if ~isreal(rawdata)
error(['There are complex values in the data. Probably a wrong' ...
' transformation'])
end
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
offset = npar-np;
fname_=M_.fname;
options_ = set_default_option(options_,'opt_gsa',1);
options_gsa_ = options_.opt_gsa;
if options_gsa_.pprior,
namfile=[fname_,'_prior'];
else
namfile=[fname_,'_mc'];
end
load([OutDir,'\',namfile],'lpmat', 'lpmat0', 'istable')
load(options_.mode_file)
%%
%%
%%
x=[lpmat0(istable,:) lpmat(istable,:)];
clear lpmat lpmat0 istable %iunstable egg yys T
B = size(x,1);
if M_.maximum_lag > 1
l1 = flipud([cumsum(M_.lead_lag_incidence(1:M_.maximum_lag-1,dr.order_var),1);ones(1,M_.endo_nbr)]);
n1 = nnz(l1);
else
n1 = M_.endo_nbr;
end
nfil=B/40;
stock_smooth = zeros(gend,n1,40);
stock_filter = zeros(gend+1,n1,40);
stock_ys = zeros(M_.endo_nbr,40);
logpo2=zeros(B,1);
%%
h = waitbar(0,'MC smoother ...');
ib=0;
ifil=0;
opt_gsa=options_.opt_gsa;
for b=1:B
ib=ib+1;
deep = x(b,:);
%deep(1:offset) = xparam1(1:offset);
logpo2(b,1) = DsgeLikelihood(deep',gend,data);
if opt_gsa.lik_only==0,
[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff] = DsgeSmoother(deep,gend,data);
stock_smooth(:,:,ib)=atT';
stock_filter(:,:,ib)=filtered_state_vector';
stock_ys(:,ib)=ys;
if ib==40,
ib=0;
ifil=ifil+1;
save([OutDir,'\',namfile,'_',num2str(ifil)],'stock_smooth','stock_filter','stock_ys')
stock_smooth = zeros(gend,n1,40);
stock_filter = zeros(gend+1,n1,40);
stock_ys = zeros(M_.endo_nbr,40);
end
end
waitbar(b/B,h,['MC smoother ...',num2str(b),'/',num2str(B)]);
end
close(h)
if opt_gsa.lik_only==0,
if ib>0,
ifil=ifil+1;
stock_smooth = stock_smooth(:,:,1:ib);
stock_filter = stock_filter(:,:,1:ib);
stock_ys = stock_ys(:,1:ib);
save([OutDir,'\',namfile,'_',num2str(ifil)],'stock_smooth','stock_filter','stock_ys')
end
end
stock_gend=gend;
stock_data=data;
save([OutDir,'\',namfile],'x','logpo2','stock_gend','stock_data','-append')

253
matlab/dynare_MC0.m
View File

@ -1,253 +0,0 @@
function dynare_MC(var_list_,OutDir)
global M_ options_ oo_ estim_params_
global bayestopt_
options_.varlist = var_list_;
options_.lgyidx2varobs = zeros(size(M_.endo_names,1),1);
for i = 1:size(M_.endo_names,1)
tmp = strmatch(deblank(M_.endo_names(i,:)),options_.varobs,'exact');
if ~isempty(tmp)
options_.lgyidx2varobs(i,1) = tmp;
end
end
options_ = set_default_option(options_,'first_obs',1);
options_ = set_default_option(options_,'prefilter',0);
options_ = set_default_option(options_,'presample',0);
options_ = set_default_option(options_,'lik_algo',1);
options_ = set_default_option(options_,'lik_init',1);
options_ = set_default_option(options_,'nograph',0);
options_ = set_default_option(options_,'mh_conf_sig',0.90);
options_ = set_default_option(options_,'mh_replic',20000);
options_ = set_default_option(options_,'mh_drop',0.5);
options_ = set_default_option(options_,'mh_jscale',0.2);
options_ = set_default_option(options_,'mh_init_scale',2*options_.mh_jscale);
options_ = set_default_option(options_,'mode_file','');
options_ = set_default_option(options_,'mode_compute',4);
options_ = set_default_option(options_,'mode_check',0);
options_ = set_default_option(options_,'prior_trunc',1e-10);
options_ = set_default_option(options_,'mh_mode',1);
options_ = set_default_option(options_,'mh_nblck',2);
options_ = set_default_option(options_,'load_mh_file',0);
options_ = set_default_option(options_,'nodiagnostic',0);
options_ = set_default_option(options_,'loglinear',0);
options_ = set_default_option(options_,'unit_root_vars',[]);
options_ = set_default_option(options_,'XTick',[]);
options_ = set_default_option(options_,'XTickLabel',[]);
options_ = set_default_option(options_,'bayesian_irf',0);
options_ = set_default_option(options_,'bayesian_th_moments',0);
options_ = set_default_option(options_,'TeX',0);
options_ = set_default_option(options_,'irf',0);
options_ = set_default_option(options_,'relative_irf',0);
options_ = set_default_option(options_,'order',1);
options_ = set_default_option(options_,'ar',5);
options_ = set_default_option(options_,'dr_algo',0);
options_ = set_default_option(options_,'linear',0);
options_ = set_default_option(options_,'drop',0);
options_ = set_default_option(options_,'replic',1);
options_ = set_default_option(options_,'hp_filter',0);
options_ = set_default_option(options_,'forecast',0);
options_ = set_default_option(options_,'smoother',0);
options_ = set_default_option(options_,'moments_varendo',0);
options_ = set_default_option(options_,'filtered_vars',0);
options_ = set_default_option(options_,'kalman_algo',1);
options_ = set_default_option(options_,'kalman_tol',10^(-12));
options_ = set_default_option(options_,'posterior_mode_estimation',1);
options_ = set_default_option(options_,'MaxNumberOfBytes',1e6);
options_ = set_default_option(options_,'xls_sheet','');
options_ = set_default_option(options_,'xls_range','');
options_ = set_default_option(options_,'steadystate_flag',0);
if exist([M_.fname '_steadystate.m'])
options_.steadystate_flag = 1;
end
%% Add something to the parser ++>
M_.dname = M_.fname; % The user should be able to choose another name
% for the directory...
pnames = [' ';'beta ';'gamm ';'norm ';'invg ';'unif ';'invg2'];
n_varobs = size(options_.varobs,1);
[xparam1,estim_params_,bayestopt_,lb,ub] = set_prior(estim_params_);
options_.mh_replic = 0;
bounds = prior_bounds(bayestopt_);
bounds(:,1)=max(bounds(:,1),lb);
bounds(:,2)=min(bounds(:,2),ub);
if any(xparam1 < bounds(:,1)) | any(xparam1 > bounds(:,2))
find(xparam1 < bounds(:,1))
find(xparam1 > bounds(:,2))
error('Initial parameter values are outside parameter bounds')
end
lb = bounds(:,1);
ub = bounds(:,2);
bayestopt_.lb = lb;
bayestopt_.ub = ub;
if ~isfield(options_,'trend_coeffs')
bayestopt_.with_trend = 0;
else
bayestopt_.with_trend = 1;
bayestopt_.trend_coeff = {};
trend_coeffs = options_.trend_coeffs;
nt = length(trend_coeffs);
for i=1:n_varobs
if i > length(trend_coeffs)
bayestopt_.trend_coeff{i} = '0';
else
bayestopt_.trend_coeff{i} = trend_coeffs{i};
end
end
end
bayestopt_.penalty = 1e8; % penalty
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
nx = nvx+nvn+ncx+ncn+np;
dr = set_state_space([]);
%% Initialization with unit-root variables
if ~isempty(options_.unit_root_vars)
n_ur = size(options_.unit_root_vars,1);
i_ur = zeros(n_ur,1);
for i=1:n_ur
i1 = strmatch(deblank(options_.unit_root_vars(i,:)),M_.endo_names(dr.order_var,:),'exact');
if isempty(i1)
error('Undeclared variable in unit_root_vars statement')
end
i_ur(i) = i1;
end
if M_.maximum_lag > 1
l1 = flipud([cumsum(M_.lead_lag_incidence(1:M_.maximum_lag-1,dr.order_var),1);ones(1,M_.endo_nbr)]);
n1 = nnz(l1);
bayestopt_.Pinf = zeros(n1,n1);
l2 = find(l1');
l3 = zeros(M_.endo_nbr,M_.maximum_lag);
l3(i_ur,:) = l1(:,i_ur)';
l3 = l3(:);
i_ur1 = find(l3(l2));
i_stable = ones(M_.endo_nbr,1);
i_stable(i_ur) = zeros(n_ur,1);
i_stable = find(i_stable);
bayestopt_.Pinf(i_ur1,i_ur1) = diag(ones(1,length(i_ur1)));
bayestopt_.i_var_stable = i_stable;
l3 = zeros(M_.endo_nbr,M_.maximum_lag);
l3(i_stable,:) = l1(:,i_stable)';
l3 = l3(:);
bayestopt_.i_T_var_stable = find(l3(l2));
else
n1 = M_.endo_nbr;
bayestopt_.Pinf = zeros(n1,n1);
bayestopt_.Pinf(i_ur,i_ur) = diag(ones(1,length(i_ur)));
l1 = ones(M_.endo_nbr,1);
l1(i_ur,:) = zeros(length(i_ur),1);
bayestopt_.i_T_var_stable = find(l1);
end
options_.lik_init = 3;
end % if ~isempty(options_.unit_root_vars)
if isempty(options_.datafile)
error('ESTIMATION: datafile option is missing')
end
if isempty(options_.varobs)
error('ESTIMATION: VAROBS is missing')
end
%% If jscale isn't specified for an estimated parameter, use
%% global option options_.jscale, set to 0.2, by default
k = find(isnan(bayestopt_.jscale));
bayestopt_.jscale(k) = options_.mh_jscale;
%% Read and demean data
rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
k = [];
k1 = [];
for i=1:n_varobs
k = [k strmatch(deblank(options_.varobs(i,:)),M_.endo_names(dr.order_var,:), ...
'exact')];
k1 = [k1 strmatch(deblank(options_.varobs(i,:)),M_.endo_names, 'exact')];
end
bayestopt_.mf = k;
bayestopt_.mfys = k1;
options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1);
gend = options_.nobs;
rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
if options_.loglinear == 1
rawdata = log(rawdata);
end
if options_.prefilter == 1
bayestopt_.mean_varobs = mean(rawdata,1);
data = transpose(rawdata-ones(gend,1)*bayestopt_.mean_varobs);
else
data = transpose(rawdata);
end
if ~isreal(rawdata)
error(['There are complex values in the data. Probably a wrong' ...
' transformation'])
end
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
offset = npar-np;
fname_=M_.fname;
options_ = set_default_option(options_,'opt_gsa',1);
options_gsa_ = options_.opt_gsa;
if options_gsa_.pprior,
namfile=[fname_,'_prior'];
else
namfile=[fname_,'_mc'];
end
load([OutDir,'\',namfile])
load(options_.mode_file)
%%
%%
%%
x=[lpmat0(istable,:) lpmat(istable,:)];
clear lpmat lpmat0 istable iunstable egg yys T
B = size(x,1);
if M_.maximum_lag > 1
l1 = flipud([cumsum(M_.lead_lag_incidence(1:M_.maximum_lag-1,dr.order_var),1);ones(1,M_.endo_nbr)]);
n1 = nnz(l1);
else
n1 = M_.endo_nbr;
end
stock_smooth = zeros(gend,n1,B);
stock_filter = zeros(gend+1,n1,B);
stock_ys = zeros(M_.endo_nbr,B);
%%
h = waitbar(0,'MC smoother ...');
for b=1:B
deep = x(b,:);
%deep(1:offset) = xparam1(1:offset);
logpo2(b,1) = DsgeLikelihood(deep',gend,data);
[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff] = DsgeSmoother(deep,gend,data);
stock_smooth(:,:,b)=atT';
stock_filter(:,:,b)=filtered_state_vector';
stock_ys(:,b)=ys;
waitbar(b/B,h,['MC smoother ...',num2str(b),'/',num2str(B)]);
end
close(h)
stock_gend=gend;
stock_data=data;
save([OutDir,'\',namfile],'x','logpo2','stock_smooth','stock_filter','stock_ys','stock_gend','stock_data','-append')

210
matlab/dynare_sensitivity.m
View File

@ -1,210 +0,0 @@
function x0=dynare_sensitivity()
% copyright Marco Ratto 2006
global M_ options_ oo_ bayestopt_
fname_ = M_.fname;
lgy_ = M_.endo_names;
x0=[];
options_ = set_default_option(options_,'opt_gsa',1);
options_gsa_ = options_.opt_gsa;
% map stability
options_gsa_ = set_default_option(options_gsa_,'stab',1);
options_gsa_ = set_default_option(options_gsa_,'redform',0);
options_gsa_ = set_default_option(options_gsa_,'pprior',1);
options_gsa_ = set_default_option(options_gsa_,'ppost',0);
options_gsa_ = set_default_option(options_gsa_,'ilptau',1);
options_gsa_ = set_default_option(options_gsa_,'morris',0);
options_gsa_ = set_default_option(options_gsa_,'morris_nliv',6);
options_gsa_ = set_default_option(options_gsa_,'morris_ntra',20);
options_gsa_ = set_default_option(options_gsa_,'Nsam',2048);
options_gsa_ = set_default_option(options_gsa_,'load_stab',0);
options_gsa_ = set_default_option(options_gsa_,'alpha2_stab',0.3);
options_gsa_ = set_default_option(options_gsa_,'ksstat',0.1);
%options_gsa_ = set_default_option(options_gsa_,'load_mh',0);
OutputDirectoryName = CheckPath('GSA');
if options_gsa_.morris,
options_gsa_.pprior=1;
options_gsa_.stab=1;
end
options_.opt_gsa = options_gsa_;
if options_gsa_.stab & ~options_gsa_.ppost,
x0 = stab_map_(options_gsa_.Nsam, options_gsa_.load_stab, options_gsa_.ksstat, options_gsa_.alpha2_stab, ...
options_gsa_.redform, options_gsa_.pprior, options_gsa_.ilptau, OutputDirectoryName);
end
% reduced form
% redform_map(namendo, namlagendo, namexo, icomp, pprior, ilog, threshold)
options_gsa_ = set_default_option(options_gsa_,'load_redform',0);
options_gsa_ = set_default_option(options_gsa_,'logtrans_redform',0);
options_gsa_ = set_default_option(options_gsa_,'threshold_redform',[]);
options_gsa_ = set_default_option(options_gsa_,'ksstat_redform',0.1);
options_gsa_ = set_default_option(options_gsa_,'alpha2_redform',0.3);
options_gsa_ = set_default_option(options_gsa_,'namendo',[]);
options_gsa_ = set_default_option(options_gsa_,'namlagendo',[]);
options_gsa_ = set_default_option(options_gsa_,'namexo',[]);
options_.opt_gsa = options_gsa_;
if options_gsa_.redform & ~isempty(options_gsa_.namendo) & ~options_gsa_.ppost,
redform_map(options_gsa_.namendo, options_gsa_.namlagendo, options_gsa_.namexo, ...
options_gsa_.load_redform, options_gsa_.pprior, options_gsa_.logtrans_redform, ...
options_gsa_.threshold_redform, options_gsa_.ksstat_redform, ...
options_gsa_.alpha2_redform, OutputDirectoryName);
end
% RMSE mapping
% function [rmse_MC, ixx] = filt_mc_(vvarvecm, loadSA, pfilt, alpha, alpha2)
options_gsa_ = set_default_option(options_gsa_,'rmse',0);
options_gsa_ = set_default_option(options_gsa_,'lik_only',0);
options_gsa_ = set_default_option(options_gsa_,'var_rmse',options_.varobs);
options_gsa_ = set_default_option(options_gsa_,'load_rmse',0);
options_gsa_ = set_default_option(options_gsa_,'pfilt_rmse',0.1);
options_gsa_ = set_default_option(options_gsa_,'istart_rmse',1);
options_gsa_ = set_default_option(options_gsa_,'alpha_rmse',0.002);
options_gsa_ = set_default_option(options_gsa_,'alpha2_rmse',1);
options_.opt_gsa = options_gsa_;
if options_gsa_.rmse,
if ~options_gsa_.ppost
if options_gsa_.pprior
a=whos('-file',[OutputDirectoryName,'\',fname_,'_prior'],'logpo2');
else
a=whos('-file',[OutputDirectoryName,'\',fname_,'_mc'],'logpo2');
end
if isempty(a),
dynare_MC([],OutputDirectoryName);
options_gsa_.load_rmse=0;
end
end
clear a;
filt_mc_(options_gsa_.var_rmse, options_gsa_.load_rmse, options_gsa_.pfilt_rmse, ...
options_gsa_.alpha_rmse, options_gsa_.alpha2_rmse, OutputDirectoryName, ...
options_gsa_.istart_rmse);
end
options_gsa_ = set_default_option(options_gsa_,'glue',0);
if options_gsa_.glue,
dr_ = oo_.dr;
if options_gsa_.ppost
load([OutputDirectoryName,'\',fname_,'_post']);
DirectoryName = CheckPath('metropolis');
else
if options_gsa_.pprior
load([OutputDirectoryName,'\',fname_,'_prior']);
else
load([OutputDirectoryName,'\',fname_,'_mc']);
end
end
nruns=size(x,1);
gend = options_.nobs;
rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
if options_.loglinear == 1
rawdata = log(rawdata);
end
if options_.prefilter == 1
%data = transpose(rawdata-ones(gend,1)*bayestopt_.mean_varobs);
data = transpose(rawdata-ones(gend,1)*mean(rawdata,1));
else
data = transpose(rawdata);
end
Obs.data = data;
Obs.time = [1:gend];
Obs.num = gend;
for j=1:size(options_.varobs,1)
Obs.name{j} = deblank(options_.varobs(j,:));
vj=deblank(options_.varobs(j,:));
jxj = strmatch(vj,lgy_(dr_.order_var,:),'exact');
js = strmatch(vj,lgy_,'exact');
if ~options_gsa_.ppost
y0=zeros(gend+1,nruns);
nb = size(stock_filter,3);
y0 = squeeze(stock_filter(:,jxj,:)) + ...
kron(stock_ys(js,:),ones(size(stock_filter,1),1));
Out(j).data = y0';
Out(j).time = [1:size(y0,1)];
else
Out(j).data = jxj;
Out(j).time = [pwd,'\',DirectoryName];
end
Out(j).name = vj;
Out(j).ini = 'yes';
Lik(j).name = ['rmse_',vj];
Lik(j).ini = 'yes';
Lik(j).isam = 1;
Lik(j).data = rmse_MC(:,j)';
if ~options_gsa_.ppost
y0 = squeeze(stock_smooth(:,jxj,:)) + ...
kron(stock_ys(js,:),ones(size(stock_smooth,1),1));
Out1(j).name = vj;
Out1(j).ini = 'yes';
Out1(j).time = [1:size(y0,1)];
Out1(j).data = y0';
else
Out1=Out;
end
ismoo(j)=jxj;
end
jsmoo = size(options_.varobs,1);
for j=1:M_.endo_nbr,
if ~ismember(j,ismoo),
jsmoo=jsmoo+1;
vj=deblank(M_.endo_names(dr_.order_var(j),:));
if ~options_gsa_.ppost
y0 = squeeze(stock_smooth(:,j,:)) + ...
kron(stock_ys(j,:),ones(size(stock_smooth,1),1));
Out1(jsmoo).time = [1:size(y0,1)];
Out1(jsmoo).data = y0';
else
Out1(jsmoo).data = j;
Out1(jsmoo).time = [pwd,'\',DirectoryName];
end
Out1(jsmoo).name = vj;
Out1(jsmoo).ini = 'yes';
end
end
tit(M_.exo_names_orig_ord,:) = M_.exo_names;
for j=1:M_.exo_nbr,
Exo(j).name = deblank(tit(j,:));
end
if ~options_gsa_.ppost
Lik(size(options_.varobs,1)+1).name = 'logpo';
Lik(size(options_.varobs,1)+1).ini = 'yes';
Lik(size(options_.varobs,1)+1).isam = 1;
Lik(size(options_.varobs,1)+1).data = -logpo2;
end
Sam.name = bayestopt_.name;
Sam.dim = [size(x) 0];
Sam.data = [x];
Rem.id = 'Original';
Rem.ind= [1:size(x,1)];
if options_gsa_.ppost
Info.dynare=M_.fname;
Info.order_var=dr_.order_var;
Out=Out1;
save([OutputDirectoryName,'\',fname_,'_post_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
%save([fname_,'_post_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info')
else
if options_gsa_.pprior
save([OutputDirectoryName,'\',fname_,'_prior_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
Out=Out1;
save([OutputDirectoryName,'\',fname_,'_prior_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
else
save([OutputDirectoryName,'\',fname_,'_mc_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
Out=Out1;
save([OutputDirectoryName,'\',fname_,'_mc_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
end
end
end

19
matlab/facDefault.m
View File

@ -1,19 +0,0 @@
function facDefault
% restores factory defined defaults
set(0,'DefaultFigurePaperType','a4')
%set(0,'DefaultFigurePaperOrientation','landscape')
set(0,'DefaultFigureColor',[0.8 0.8 0.8])
set(0,'DefaultTextColor',[0 0 0])
set(0,'DefaultAxesUnits','normalized')
set(0,'DefaultAxesLineStyleOrder', ...
{'-',':','--','-.'})
set(0,'DefaultAxesColor',[1 1 1],...
'DefaultAxesXColor', [0 0 0], ...
'DefaultAxesYColor', [0 0 0], ...
'DefaultAxesZColor', [0 0 0])
set(0, ...
'DefaultAxesColorOrder', ...
[0 0 0; 0 0 1; 0 1 0; 1 0 0; 1 0 1; 0 1 1; 1 1 0])
set(0,'DefaultPatchEdgeColor','factory')

643
matlab/filt_mc_.m
View File

@ -1,643 +0,0 @@
function [rmse_MC, ixx] = filt_mc_(vvarvecm, loadSA, pfilt, alpha, alpha2, OutDir, istart, alphaPC)
% copyright Marco Ratto 2006
global bayestopt_ estim_params_ M_ options_ oo_
if nargin<1 | isempty(vvarvecm),
vvarvecm = options_.varobs;
end
if nargin<2,
loadSA=0;
end
if nargin<3 | isempty(pfilt),
pfilt=0.1; % cut the best 10% of runs
end
if nargin<4 | isempty(alpha),
alpha=0.002;
end
if nargin<5 | isempty(alpha2),
alpha2=0.5;
end
if nargin<7 | isempty(istart),
istart=1;
end
if nargin<8,
alphaPC=0.5;
end
fname_ = M_.fname;
lgy_ = M_.endo_names;
dr_ = oo_.dr;
disp(' ')
disp(' ')
disp('Starting sensitivity analysis')
disp('for the fit of EACH observed series ...')
disp(' ')
disp('Deleting old SA figures...')
a=dir([OutDir,'\*.*']);
if options_.opt_gsa.ppost,
tmp=['_SA_fit_post'];
else
if options_.opt_gsa.pprior
tmp=['_SA_fit_prior'];
else
tmp=['_SA_fit_mc'];
end
end
for j=1:length(a),
if strmatch([fname_,tmp],a(j).name),
disp(a(j).name)
delete([OutDir,'\',a(j).name])
end,
end
disp('done !')
nshock=estim_params_.nvx + estim_params_.nvn + estim_params_.ncx + estim_params_.ncn;
npar=estim_params_.np;
for j=1:npar+nshock,
if j>nshock
if isfield(oo_,'posterior_mode'),
xparam1(j)=oo_.posterior_mode.parameters.(bayestopt_.name{j});
end
if isfield(oo_,'posterior_mean'),
xparam1_mean(j)=oo_.posterior_mean.parameters.(bayestopt_.name{j});
end
else
if isfield(oo_,'posterior_mode'),
xparam1(j)=oo_.posterior_mode.shocks_std.(bayestopt_.name{j});
end
if isfield(oo_,'posterior_mean'),
xparam1_mean(j)=oo_.posterior_mean.shocks_std.(bayestopt_.name{j});
end
end
end
if options_.opt_gsa.ppost,
fnamtmp=[fname_,'_post'];
DirectoryName = CheckPath('metropolis');
else
if options_.opt_gsa.pprior
fnamtmp=[fname_,'_prior'];
else
fnamtmp=[fname_,'_mc'];
end
end
if ~loadSA,
if exist('xparam1','var')
set_all_parameters(xparam1);
steady_;
ys_mode=oo_.steady_state;
end
if exist('xparam1_mean','var')
set_all_parameters(xparam1_mean);
steady_;
ys_mean=oo_.steady_state;
end
eval(options_.datafile)
if ~options_.opt_gsa.ppost
load([OutDir,'\',fnamtmp]);
else
load([DirectoryName '/' M_.fname '_data.mat']);
filfilt = dir([DirectoryName '/' M_.fname '_filter*.mat']);
filparam = dir([DirectoryName '/' M_.fname '_param*.mat']);
x=[];
logpo2=[];
sto_ys=[];
for j=1:length(filparam),
%load([DirectoryName '/' M_.fname '_param',int2str(j),'.mat']);
if isempty(strmatch([M_.fname '_param_irf'],filparam(j).name))
load([DirectoryName '/' filparam(j).name]);
x=[x; stock];
logpo2=[logpo2; stock_logpo];
sto_ys=[sto_ys; stock_ys];
clear stock stock_logpo stock_ys;
end
end
logpo2=-logpo2;
end
nruns=size(x,1);
nfilt=floor(pfilt*nruns);
if options_.opt_gsa.ppost | (options_.opt_gsa.ppost==0 & options_.opt_gsa.lik_only==0)
disp(' ')
disp('Computing RMSE''s...')
fobs = options_.first_obs;
nobs=options_.nobs;
for i=1:size(vvarvecm,1),
vj=deblank(vvarvecm(i,:));
if options_.prefilter == 1
%eval([vj,'=',vj,'-bayestopt_.mean_varobs(i);'])
eval([vj,'=',vj,'-mean(',vj,',1);'])
end
jxj = strmatch(vj,lgy_(dr_.order_var,:),'exact');
js = strmatch(vj,lgy_,'exact');
if exist('xparam1','var')
eval(['rmse_mode(i) = sqrt(mean((',vj,'(fobs-1+istart:fobs-1+nobs)-oo_.steady_state(js)-oo_.FilteredVariables.',vj,'(istart:end-1)).^2));'])
end
y0=zeros(nobs+1,nruns);
if options_.opt_gsa.ppost
%y0=zeros(nobs+max(options_.filter_step_ahead),nruns);
nbb=0;
for j=1:length(filfilt),
load([DirectoryName '/' M_.fname '_filter',num2str(j),'.mat']);
nb = size(stock,4);
% y0(:,nbb+1:nbb+nb)=squeeze(stock(1,jxj,:,:)) + ...
% kron(sto_ys(nbb+1:nbb+nb,js)',ones(size(stock,3),1));
y0(:,nbb+1:nbb+nb)=squeeze(stock(1,jxj,1:nobs+1,:)) + ...
kron(sto_ys(nbb+1:nbb+nb,js)',ones(nobs+1,1));
%y0(:,:,size(y0,3):size(y0,3)+size(stock,3))=stock;
nbb=nbb+nb;
clear stock;
end
else
y0 = squeeze(stock_filter(:,jxj,:)) + ...
kron(stock_ys(js,:),ones(size(stock_filter,1),1));
end
y0M=mean(y0,2);
for j=1:nruns,
eval(['rmse_MC(j,i) = sqrt(mean((',vj,'(fobs-1+istart:fobs-1+nobs)-y0(istart:end-1,j)).^2));'])
end
if exist('xparam1_mean','var')
%eval(['rmse_pmean(i) = sqrt(mean((',vj,'(fobs-1+istart:fobs-1+nobs)-y0M(istart:end-1)).^2));'])
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1_mean,stock_gend,stock_data);
y0 = ahat(jxj,:)' + ...
kron(ys_mean(js,:),ones(size(ahat,2),1));
eval(['rmse_pmean(i) = sqrt(mean((',vj,'(fobs-1+istart:fobs-1+nobs)-y0(istart:end-1)).^2));'])
end
end
clear stock_filter;
end
for j=1:nruns,
lnprior(j,1) = priordens(x(j,:),bayestopt_.pshape,bayestopt_.p1,bayestopt_.p2,bayestopt_.p3,bayestopt_.p4);
end
likelihood=logpo2(:)+lnprior(:);
disp('... done!')
if options_.opt_gsa.ppost
save([OutDir,'\',fnamtmp], 'x', 'logpo2', 'likelihood', 'rmse_MC', 'rmse_mode','rmse_pmean')
else
if options_.opt_gsa.lik_only
save([OutDir,'\',fnamtmp], 'likelihood', '-append')
else
save([OutDir,'\',fnamtmp], 'likelihood', 'rmse_MC', 'rmse_mode','rmse_pmean','-append')
end
end
else
if options_.opt_gsa.lik_only & options_.opt_gsa.ppost==0
load([OutDir,'\',fnamtmp],'x','logpo2','likelihood');
else
load([OutDir,'\',fnamtmp],'x','logpo2','likelihood','rmse_MC','rmse_mode','rmse_pmean');
end
lnprior=likelihood(:)-logpo2(:);
nruns=size(x,1);
nfilt=floor(pfilt*nruns);
end
% smirnov tests
nfilt0=nfilt*ones(size(vvarvecm,1),1);
logpo2=logpo2(:);
if ~options_.opt_gsa.ppost
[dum, ipost]=sort(logpo2);
[dum, ilik]=sort(likelihood);
end
if ~options_.opt_gsa.ppost & options_.opt_gsa.lik_only
if options_.opt_gsa.pprior
anam='SA_fit_prior_post';
else
anam='SA_fit_mc_post';
end
stab_map_1(x, ipost(1:nfilt), ipost(nfilt+1:end), anam, 1,[],OutDir);
stab_map_2(x(ipost(1:nfilt),:),alpha2,anam, OutDir);
if options_.opt_gsa.pprior
anam='SA_fit_prior_lik';
else
anam='SA_fit_mc_lik';
end
stab_map_1(x, ilik(1:nfilt), ilik(nfilt+1:end), anam, 1,[],OutDir);
stab_map_2(x(ilik(1:nfilt),:),alpha2,anam, OutDir);
else
for i=1:size(vvarvecm,1),
[dum, ixx(:,i)]=sort(rmse_MC(:,i));
if options_.opt_gsa.ppost,
%nfilt0(i)=length(find(rmse_MC(:,i)<rmse_pmean(i)));
rmse_txt=rmse_pmean;
else
if options_.opt_gsa.pprior,
rmse_txt=rmse_mode;
else
%nfilt0(i)=length(find(rmse_MC(:,i)<rmse_pmean(i)));
rmse_txt=rmse_pmean;
end
end
for j=1:npar+nshock,
[H,P,KSSTAT] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j), alpha);
[H1,P1,KSSTAT1] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,1);
[H2,P2,KSSTAT2] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,-1);
if H1 & H2==0,
SS(j,i)=1;
elseif H1==0,
SS(j,i)=-1;
else
SS(j,i)=0;
end
PP(j,i)=P;
end
end
ifig=0;
for i=1:size(vvarvecm,1),
if mod(i,9)==1,
ifig=ifig+1;
figure('name',['Prior ',int2str(ifig)])
end
subplot(3,3,i-9*(ifig-1))
h=cumplot(lnprior(ixx(1:nfilt0(i),i)));
set(h,'color','red')
hold on, cumplot(lnprior)
h=cumplot(lnprior(ixx(nfilt0(i)+1:end,i)));
set(h,'color','green')
title(vvarvecm(i,:))
if mod(i,9)==0 | i==size(vvarvecm,1)
if options_.opt_gsa.ppost
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_post_lnprior',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_post_lnprior',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_post_lnprior',int2str(ifig)]);
else
if options_.opt_gsa.pprior
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_prior_lnprior',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_prior_lnprior',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_prior_lnprior',int2str(ifig)]);
else
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_mc_lnprior',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_mc_lnprior',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_mc_lnprior',int2str(ifig)]);
end
end
close(gcf)
end
end
ifig=0;
for i=1:size(vvarvecm,1),
if mod(i,9)==1,
ifig=ifig+1;
figure('name',['Likelihood ',int2str(ifig)])
end
subplot(3,3,i-9*(ifig-1))
h=cumplot(likelihood(ixx(1:nfilt0(i),i)));
set(h,'color','red')
hold on, h=cumplot(likelihood);
h=cumplot(likelihood(ixx(nfilt0(i)+1:end,i)));
set(h,'color','green')
title(vvarvecm(i,:))
if mod(i,9)==0 | i==size(vvarvecm,1)
if options_.opt_gsa.ppost
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_post_lnlik',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_post_lnlik',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_post_lnlik',int2str(ifig)]);
else
if options_.opt_gsa.pprior
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_prior_lnlik',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_prior_lnlik',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_prior_lnlik',int2str(ifig)]);
else
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_mc_lnlik',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_mc_lnlik',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_mc_lnlik',int2str(ifig)]);
end
end
close(gcf)
end
end
ifig=0;
for i=1:size(vvarvecm,1),
if mod(i,9)==1,
ifig=ifig+1;
figure('name',['Posterior ',int2str(ifig)])
end
subplot(3,3,i-9*(ifig-1))
h=cumplot(logpo2(ixx(1:nfilt0(i),i)));
set(h,'color','red')
hold on, h=cumplot(logpo2);
h=cumplot(logpo2(ixx(nfilt0(i)+1:end,i)));
set(h,'color','green')
title(vvarvecm(i,:))
if mod(i,9)==0 | i==size(vvarvecm,1)
if options_.opt_gsa.ppost
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_post_lnpost',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_post_lnpost',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_post_lnpost',int2str(ifig)]);
else
if options_.opt_gsa.pprior
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_prior_lnpost',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_prior_lnpost',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_prior_lnpost',int2str(ifig)]);
else
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_mc_lnpost',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_mc_lnpost',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_mc_lnpost',int2str(ifig)]);
end
end
close(gcf)
end
end
param_names='';
for j=1:npar+nshock,
param_names=str2mat(param_names, bayestopt_.name{j});
end
param_names=param_names(2:end,:);
disp(' ')
disp('RMSE over the MC sample:')
disp(' min yr RMSE max yr RMSE')
for j=1:size(vvarvecm,1),
disp([vvarvecm(j,:), sprintf('%15.5g',[(min(rmse_MC(:,j))) [(max(rmse_MC(:,j)))]])])
end
invar = find( std(rmse_MC)./mean(rmse_MC)<=0.0001 );
if ~isempty(invar)
disp(' ')
disp(' ')
disp('RMSE is not varying significantly over the MC sample for the following variables:')
disp(vvarvecm(invar,:))
disp('These variables are excluded from SA')
disp('[Unless you treat these series as exogenous, there is something wrong in your estimation !]')
end
ivar = find( std(rmse_MC)./mean(rmse_MC)>0.0001 );
vvarvecm=vvarvecm(ivar,:);
rmse_MC=rmse_MC(:,ivar);
disp(' ')
if options_.opt_gsa.ppost==0 & options_.opt_gsa.pprior,
disp(['Sample filtered the ',num2str(pfilt*100),'% best RMSE''s for each observed series ...' ])
else
disp(['Sample filtered the best RMSE''s smaller than RMSE at the posterior mean ...' ])
end
% figure, boxplot(rmse_MC)
% set(gca,'xticklabel',vvarvecm)
% saveas(gcf,[fname_,'_SA_RMSE'])
disp(' ')
disp(' ')
disp('RMSE ranges after filtering:')
if options_.opt_gsa.ppost==0 & options_.opt_gsa.pprior,
disp([' best ',num2str(pfilt*100),'% filtered remaining 90%'])
disp([' min max min max posterior mode'])
else
disp([' best filtered remaining '])
disp([' min max min max posterior mean'])
end
for j=1:size(vvarvecm,1),
disp([vvarvecm(j,:), sprintf('%15.5g',[min(rmse_MC(ixx(1:nfilt0(j),j),j)) ...
max(rmse_MC(ixx(1:nfilt0(j),j),j)) ...
min(rmse_MC(ixx(nfilt0(j)+1:end,j),j)) ...
max(rmse_MC(ixx(nfilt0(j)+1:end,j),j)) ...
rmse_txt(j)])])
% disp([vvarvecm(j,:), sprintf('%15.5g',[min(logpo2(ixx(1:nfilt,j))) ...
% max(logpo2(ixx(1:nfilt,j))) ...
% min(logpo2(ixx(nfilt+1:end,j))) ...
% max(logpo2(ixx(nfilt+1:end,j)))])])
end
SP=zeros(npar+nshock,size(vvarvecm,1));
for j=1:size(vvarvecm,1),
ns=find(PP(:,j)<alpha);
SP(ns,j)=ones(size(ns));
SS(:,j)=SS(:,j).*SP(:,j);
end
for j=1:npar+nshock, %estim_params_.np,
nsp(j)=length(find(SP(j,:)));
end
snam0=param_names(find(nsp==0),:);
snam1=param_names(find(nsp==1),:);
snam2=param_names(find(nsp>1),:);
snam=param_names(find(nsp>0),:);
% snam0=bayestopt_.name(find(nsp==0));
% snam1=bayestopt_.name(find(nsp==1));
% snam2=bayestopt_.name(find(nsp>1));
% snam=bayestopt_.name(find(nsp>0));
nsnam=(find(nsp>1));
disp(' ')
disp(' ')
disp('These parameters do not affect significantly the fit of ANY observed series:')
disp(snam0)
disp(' ')
disp('These parameters affect ONE single observed series:')
disp(snam1)
disp(' ')
disp('These parameters affect MORE THAN ONE observed series: trade off exists!')
disp(snam2)
%pnam=bayestopt_.name(end-estim_params_.np+1:end);
pnam=bayestopt_.name;
% plot trade-offs
a00=jet(size(vvarvecm,1));
for ix=1:ceil(length(nsnam)/6),
figure,
for j=1+6*(ix-1):min(size(snam2,1),6*ix),
subplot(2,3,j-6*(ix-1))
%h0=cumplot(x(:,nsnam(j)+nshock));
h0=cumplot(x(:,nsnam(j)));
set(h0,'color',[0 0 0])
hold on,
np=find(SP(nsnam(j),:));
%a0=jet(nsp(nsnam(j)));
a0=a00(np,:);
for i=1:nsp(nsnam(j)), %size(vvarvecm,1),
%h0=cumplot(x(ixx(1:nfilt,np(i)),nsnam(j)+nshock));
h0=cumplot(x(ixx(1:nfilt0(np(i)),np(i)),nsnam(j)));
set(h0,'color',a0(i,:))
end
ydum=get(gca,'ylim');
%xdum=xparam1(nshock+nsnam(j));
xdum=xparam1(nsnam(j));
h1=plot([xdum xdum],ydum);
set(h1,'color',[0.85 0.85 0.85],'linewidth',2)
h0=legend(str2mat('base',vvarvecm(np,:)),0);
set(h0,'fontsize',6)
%h0=legend({'base',vnam{np}}',0);
xlabel('')
set(findobj(get(h0,'children'),'type','text'),'interpreter','none')
title([pnam{nsnam(j)}],'interpreter','none')
if options_.opt_gsa.ppost
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_post_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_post_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_post_' int2str(ix)]);
else
if options_.opt_gsa.pprior
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_prior_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_prior_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_prior_' int2str(ix)]);
else
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_mc_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_mc_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_mc_' int2str(ix)]);
end
end
end
end
close all
for j=1:size(SP,2),
nsx(j)=length(find(SP(:,j)));
end
number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two.
bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton.
%kernel_function = 'uniform'; % Gaussian kernel for Fast Fourrier Transform approximaton.
for ix=1:ceil(length(nsnam)/6),
figure,
for j=1+6*(ix-1):min(size(snam2,1),6*ix),
subplot(2,3,j-6*(ix-1))
optimal_bandwidth = mh_optimal_bandwidth(x(:,nsnam(j)),size(x,1),bandwidth,kernel_function);
[x1,f1] = kernel_density_estimate(x(:,nsnam(j)),number_of_grid_points,...
optimal_bandwidth,kernel_function);
h0 = plot(x1, f1,'k');
hold on,
np=find(SP(nsnam(j),:));
%a0=jet(nsp(nsnam(j)));
a0=a00(np,:);
for i=1:nsp(nsnam(j)), %size(vvarvecm,1),
optimal_bandwidth = mh_optimal_bandwidth(x(ixx(1:nfilt0(np(i)),np(i)),nsnam(j)),nfilt,bandwidth,kernel_function);
[x1,f1] = kernel_density_estimate(x(ixx(1:nfilt0(np(i)),np(i)),nsnam(j)),number_of_grid_points,...
optimal_bandwidth,kernel_function);
h0 = plot(x1, f1);
set(h0,'color',a0(i,:))
end
ydum=get(gca,'ylim');
set(gca,'ylim',[0 ydum(2)]);
%xdum=xparam1(nshock+nsnam(j));
xdum=xparam1(nsnam(j));
h1=plot([xdum xdum],[0 ydum(2)]);
set(h1,'color',[0.85 0.85 0.85],'linewidth',2)
h0=legend(str2mat('base',vvarvecm(np,:)),0);
set(h0,'fontsize',6)
%h0=legend({'base',vnam{np}}',0);
xlabel('')
set(findobj(get(h0,'children'),'type','text'),'interpreter','none')
title([pnam{nsnam(j)}],'interpreter','none')
if options_.opt_gsa.ppost
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_post_dens_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_post_dens_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_post_dens_' int2str(ix)]);
else
if options_.opt_gsa.pprior
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_prior_dens_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_prior_dens_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_prior_dens_' int2str(ix)]);
else
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_mc_dens_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_mc_dens_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_mc_dens_' int2str(ix)]);
end
end
end
end
close all
% for j=1:size(SP,2),
% nfig=0;
% np=find(SP(:,j));
% for i=1:nsx(j), %size(vvarvecm,1),
% if mod(i,12)==1,
% nfig=nfig+1;
% %figure('name',['Sensitivity of fit of ',vnam{j}]),
% figure('name',['Sensitivity of fit of ',deblank(vvarvecm(j,:)),' ',num2str(nfig)]),
% end
%
% subplot(3,4,i-12*(nfig-1))
% optimal_bandwidth = mh_optimal_bandwidth(x(ixx(1:nfilt,j),np(i)),nfilt,bandwidth,kernel_function);
% [x1,f1] = kernel_density_estimate(x(ixx(1:nfilt,j),np(i)),number_of_grid_points,...
% optimal_bandwidth,kernel_function);
% plot(x1, f1,':k','linewidth',2)
% optimal_bandwidth = mh_optimal_bandwidth(x(ixx(nfilt+1:end,j),np(i)),nruns-nfilt,bandwidth,kernel_function);
% [x1,f1] = kernel_density_estimate(x(ixx(nfilt+1:end,j),np(i)),number_of_grid_points,...
% optimal_bandwidth,kernel_function);
% hold on, plot(x1, f1,'k','linewidth',2)
% ydum=get(gca,'ylim');
% %xdum=xparam1(nshock+np(i));
% xdum=xparam1(np(i));
% h1=plot([xdum xdum],ydum);
% set(h1,'color',[0.85 0.85 0.85],'linewidth',2)
% %xdum1=mean(x(ixx(1:nfilt,j),np(i)+nshock));
% xdum1=mean(x(ixx(1:nfilt,j),np(i)));
% h2=plot([xdum1 xdum1],ydum);
% set(h2,'color',[0 1 0],'linewidth',2)
% % h0=cumplot(x(nfilt+1:end,np(i)+nshock));
% % set(h0,'color',[1 1 1])
% % hold on,
% % h0=cumplot(x(ixx(1:nfilt,j),np(i)+nshock));
% % set(h0,'linestyle',':','color',[1 1 1])
% %title([pnam{np(i)}])
% title([pnam{np(i)},'. K-S prob ', num2str(PP(np(i),j))],'interpreter','none')
% xlabel('')
% if mod(i,12)==0 | i==nsx(j),
% saveas(gcf,[fname_,'_SA_fit_',deblank(vvarvecm(j,:)),'_',int2str(nfig)])
% close(gcf)
% end
% end
% end
disp(' ')
disp(' ')
disp('Sensitivity table (significance and direction):')
vav=char(zeros(1, size(param_names,2)+3 ));
ibl = 12-size(vvarvecm,2);
for j=1:size(vvarvecm,1),
vav = [vav, char(zeros(1,ibl)),vvarvecm(j,:)];
end
disp(vav)
for j=1:npar+nshock, %estim_params_.np,
%disp([param_names(j,:), sprintf('%8.5g',SP(j,:))])
disp([param_names(j,:),' ', sprintf('%12.3g',PP(j,:))])
disp([char(zeros(1, size(param_names,2)+3 )),sprintf(' (%6g)',SS(j,:))])
end
disp(' ')
disp(' ')
disp('Starting bivariate analysis:')
for i=1:size(vvarvecm,1)
if options_.opt_gsa.ppost
fnam = ['SA_fit_post_',deblank(vvarvecm(i,:))];
else
if options_.opt_gsa.pprior
fnam = ['SA_fit_prior_',deblank(vvarvecm(i,:))];
else
fnam = ['SA_fit_mc_',deblank(vvarvecm(i,:))];
end
end
stab_map_2(x(ixx(1:nfilt0(i),i),:),alpha2,fnam, OutDir);
% [pc,latent,explained] = pcacov(c0);
% %figure, bar([explained cumsum(explained)])
% ifig=0;
% j2=0;
% for j=1:npar+nshock,
% i2=find(abs(pc(:,j))>alphaPC);
% if ~isempty(i2),
% j2=j2+1;
% if mod(j2,12)==1,
% ifig=ifig+1;
% figure('name',['PCA of the filtered sample ',deblank(vvarvecm(i,:)),' ',num2str(ifig)]),
% end
% subplot(3,4,j2-(ifig-1)*12)
% bar(pc(i2,j)),
% set(gca,'xticklabel',bayestopt_.name(i2)),
% set(gca,'xtick',[1:length(i2)])
% title(['PC ',num2str(j),'. Explained ',num2str(explained(j)),'%'])
% end
% if (mod(j2,12)==0 | j==(npar+nshock)) & j2,
% saveas(gcf,[fname_,'_SA_PCA_',deblank(vvarvecm(i,:)),'_',int2str(ifig)])
% end
% end
% close all
end
end

48
matlab/fiscal_stab.m
View File

@ -1,48 +0,0 @@
function fiscal_stab(th_std_mh, vname, tname, aname)
global M_
clear fisc
if nargin<4, aname=''; end
fname_=M_.fname;
nv = size(th_std_mh,1);
for j=2:size(th_std_mh,3),
figure
fisc=squeeze((th_std_mh(:,:,j)-th_std_mh(:,:,1))./th_std_mh(:,:,1));
fisc=fisc';
fiscpol(:,j-1)=mean(fisc)';
fiscpol_std(:,j-1)=std(fisc)';
fisc=fisc.*100;
fy = fisc;
st=std(fy);
ij=find(st~=0);
%boxplot(fisc,'label',{' ',' ',' ',' ',' ',' '})
boxplot(fy(:,ij),1,'.',[],3)
hold on,
plot([1 length(ij)], [0 0], 'g')
yl=get(gca,'ylim');
for i=1:length(ij),
%text(i,(-1)^i*max((-1)^i*fisc(:,i)),cname{i})
if abs(min(fisc(:,ij(i)))-yl(1))>abs(max(fisc(:,ij(i)))-yl(2)),
text(i,min(fisc(:,ij(i))),[vname{ij(i)},' '],'rotation',90,'verticalalignment','middle','horizontalalignment','right')
else
text(i,max(fisc(:,ij(i))),[' ',vname{ij(i)}],'rotation',90,'verticalalignment','middle','horizontalalignment','left')
end
end
title(tname{j-1})
xlabel('')
ylabel('% stabilisation effect')
set(gca,'xticklabel','')
saveas(gcf,[fname_,'_fiscal',aname,num2str(j-1)])
eval(['print -depsc2 ' fname_,'_fiscal',aname,num2str(j-1)]);
eval(['print -dpdf ' fname_,'_fiscal',aname,num2str(j-1)]);
disp(tname{j-1})
sprintf('%5.2f%% (%5.2f%%) \n',[fiscpol(:,j-1)'; fiscpol_std(:,j-1)'].*100)
end
aa=squeeze(mean(th_std_mh,2));
save([fname_,'_fiscal_mh',aname], 'aa', 'fiscpol', 'fiscpol_std', '-append')

141
matlab/fstop.m
View File

@ -1,141 +0,0 @@
function [x, options] = fstop(funfcn,x,options,grad,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10)
% FSTOP Optimisation with user intervention
if nargin<3, options = []; end
options = foptions(options);
prnt = options(1);
tol = options(2);
tol2 = options(3);
evalstr = [funfcn];
stop_panel = P10;
if ~any(funfcn<48)
evalstr=[evalstr, '(x'];
for i=1:nargin - 4
evalstr = [evalstr,',P',int2str(i)];
end
evalstr = [evalstr, ')'];
end
n = prod(size(x));
if (~options(14))
options(14) = 200*n;
end
xin = x(:);
v = xin;
x(:) = v; fv = eval(evalstr);
usual_delta = 0.05;
zero_term_delta = 0.00025;
for j = 1:n
y = xin;
if y(j) ~= 0
y(j) = (1 + usual_delta)*y(j);
else
y(j) = zero_term_delta;
end
v = [v y];
x(:) = y; f = eval(evalstr);
fv = [fv f];
end
[fv,j] = sort(fv);
v = v(:,j);
cnt = n+1;
if prnt
clc
format compact
format short e
home
cnt
disp('initial ')
disp(' ')
v
f
end
alpha = 1; beta = 1/2; gamma = 2;
[n,np1] = size(v);
onesn = ones(1,n);
ot = 2:n+1;
on = 1:n;
while cnt < options(14)
if max(max(abs(v(:,ot)-v(:,onesn)))) <= tol & ...
max(abs(fv(1)-fv(ot))) <= tol2
break
end
vbar = (sum(v(:,on)')/n)';
vr = (1 + alpha)*vbar - alpha*v(:,n+1);
x(:) = vr;
fr = eval(evalstr);
cnt = cnt + 1;
vk = vr; fk = fr; how = 'reflect ';
if fr < fv(n)
if fr < fv(1)
ve = gamma*vr + (1-gamma)*vbar;
x(:) = ve;
fe = eval(evalstr);
cnt = cnt + 1;
if fe < fv(1)
vk = ve; fk = fe;
how = 'expand ';
end
end
else
vt = v(:,n+1); ft = fv(n+1);
if fr < ft
vt = vr; ft = fr;
end
vc = beta*vt + (1-beta)*vbar;
x(:) = vc;
fc = eval(evalstr);
cnt = cnt + 1;
if fc < fv(n)
vk = vc; fk = fc;
how = 'contract';
else
for j = 2:n
v(:,j) = (v(:,1) + v(:,j))/2;
x(:) = v(:,j);
fv(j) = eval(evalstr);
end
cnt = cnt + n-1;
vk = (v(:,1) + v(:,n+1))/2;
x(:) = vk;
fk = eval(evalstr);
cnt = cnt + 1;
how = 'shrink ';
end
end
v(:,n+1) = vk;
fv(n+1) = fk;
[fv,j] = sort(fv);
v = v(:,j);
if prnt
home
cnt
disp(how)
disp(' ')
v
fv
end
if ~isempty(stop_panel)
stop_panel=get(gcf,'UserData');
if stop_panel==1;
options(14)=cnt;
disp(' ')
disp(['Optimisation terminated by user after ', ...
int2str(options(14)),' function evaluations.']);
end
end
end
x(:) = v(:,1);
if prnt, format, end
options(10)=cnt;
options(8)=min(fv);
if cnt==options(14)
if options(1) >= 0
if get(gcf,'UserData')~=1
disp(['Warning: Maximum number of iterations (', ...
int2str(options(14)),') has been exceeded']);
disp( ' (increase OPTIONS(14)).')
end
end
end
% end of m-file

BIN
matlab/fstop.p
View File

Binary file not shown.

36
matlab/getFilterRuns.m
View File

@ -1,36 +0,0 @@
function y0 = getFilterRuns(varlist)
global M_ oo_
y0=[];
dirname = CheckPath('Metropolis');
filfilt = dir([dirname '/' M_.fname '_filter*.mat']);
order_var=oo_.dr.order_var;
for j=1:size(varlist,1)
jxj(j)=strmatch(deblank(varlist(j,:)),M_.endo_names,'exact');
jfilt(j)=find(order_var==jxj(j));
end
filparam = dir([dirname '/' M_.fname '_param*.mat']);
sto_ys=[];
for j=1:length(filparam),
%load([DirectoryName '/' M_.fname '_param',int2str(j),'.mat']);
if isempty(strmatch([M_.fname '_param_irf'],filparam(j).name))
load([dirname '/' filparam(j).name],'stock_ys');
sto_ys=[sto_ys; stock_ys];
clear stock stock_logpo stock_ys;
end
end
nbb=0;
for j=1:length(filfilt),
load([dirname '/' M_.fname '_filter',num2str(j),'.mat']);
nb = size(stock,4);
%y0(:,nbb+1:nbb+nb)=squeeze(stock(1,jxj,:,:));
for i=1:length(jxj)
y0(i,:,nbb+1:nbb+nb)=squeeze(stock(1,jfilt(i),:,:)) + ...
kron(sto_ys(nbb+1:nbb+nb,jxj(i))',ones(size(stock,3),1));
end
%y0(:,:,size(y0,3):size(y0,3)+size(stock,3))=stock;
nbb=nbb+nb;
end

19
matlab/getIRFRuns.m
View File

@ -1,19 +0,0 @@
function y0 = getIRFRuns(jxj,jexo,type)
global M_
if nargin<3, type='metropolis'; end,
if strcmpi(type,'posterior')
dirname = [CheckPath('metropolis') '/' ];
elseif strcmpi(type,'gsa')
dirname = [CheckPath('GSA') '/' ];
else
dirname = [CheckPath('prior') '/' ];
end
filIRF = dir([dirname '/' M_.fname '_IRFs*.mat']);
y0=[];
for file = 1:length(filIRF),
load([dirname '/' M_.fname '_IRFs' int2str(file)]);
y0 = [y0; squeeze(STOCK_IRF(:,jxj,jexo,:))];
end

15
matlab/get_mean.m
View File

@ -1,15 +0,0 @@
function y0 = get_mean(varargin);
global oo_ M_
ys_ = oo_.steady_state;
lgy_ = M_.endo_names;
lgobs_= [];
mfys=[];
for j=1:length(varargin),
dum = strmatch(varargin{j},lgy_,'exact');
mfys = [mfys dum];
end
y0 = ys_(mfys);

17
matlab/get_posterior_mode.m
View File

@ -1,17 +0,0 @@
function x=get_posterior_mode
global oo_
fn=fieldnames(oo_.posterior_density.shocks_std);
for j=1:size(fn,1),
xx=getfield(oo_.posterior_density.shocks_std,fn{j});
[dum, i]=max(xx(:,2));
x(j,1)=xx(i,1);
end
offset=length(x);
fn=fieldnames(oo_.posterior_density);
for j=1:size(fn,1)-1,
xx=getfield(oo_.posterior_density,fn{j});
[dum, i]=max(xx(:,2));
x(j+offset,1)=xx(i,1);
end

34
matlab/gr2perc.m
View File

@ -1,34 +0,0 @@
function gr2perc(cname, lname, varargin)
global M_ options_
fname_ =M_.fname;
for i=1:length(varargin)
nfig=0;
nplo=0;
for j=1:length(cname),
if evalin('base',['exist(''',cname{j},'_',varargin{i},''')']),
nplo=nplo+1;
ytemp_ = evalin('base',[cname{j},'_',varargin{i}]);
%y=exp(cumsum(ytemp_))-1;
y=cumsum(ytemp_);
if mod(nplo,9)==1,
figure('name',['Orthogonalised shocks to ',varargin{i}, ' in perc''s'])
nfig=nfig+1;
end
subplot(3,3,nplo-9*(nfig-1))
plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
hold on,
plot([1:options_.irf], y,'k'),
assignin('base',[lname{j},'_',varargin{i}],y);
title(lname{j},'interpreter','none')
hold off,
end
if mod(nplo,9)==0 | j==length(cname),
saveas(gcf,[fname_,'_IRF_',varargin{i},'_P',int2str(nfig)])
eval(['print -depsc2 ' fname_,'_IRF_',varargin{i},'_P',int2str(nfig)]);
eval(['print -dpdf ' fname_,'_IRF_',varargin{i},'_P',int2str(nfig)]);
close(gcf)
end
end
end

66
matlab/gr2perc_posterior.m
View File

@ -1,66 +0,0 @@
function gr2perc_posterior(cname, lname, varargin)
global M_ oo_ options_
fname_ =M_.fname;
DirectoryName = CheckPath('Output');
iexo=[];
for i=1:length(varargin)
iexo=[iexo, strmatch(varargin{i},M_.exo_names,'exact')];
end
iendo=[];
for j=1:length(cname),
iendo=[iendo, strmatch(cname{j},M_.endo_names,'exact')];
end
nsubplo=9;
nn=3;
y=getIRFRuns(iendo,iexo);
for i=1:length(varargin)
nfig=0;
nplo=0;
for j=1:length(cname),
y0=squeeze(y(:,j,i,:));
y1=cumsum(y0);
for k=1:options_.irf,
[MeanIRF(k,j,i),MedianIRF(k,j,i),VarIRF(k,j,i),HPDIRF(k,:,j,i),DistribIRF(k,:,j,i)] = ...
posterior_moments(y1(k,:),0);
end
name = [lname{j} '_' varargin{i}];
eval(['oo_.PosteriorIRF.Mean.' name ' = MeanIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.Median.' name ' = MedianIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.Var.' name ' = VarIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.Distribution.' name ' = DistribIRF(:,:,j,i);']);
eval(['oo_.PosteriorIRF.HPDinf.' name ' = HPDIRF(:,1,j,i);']);
eval(['oo_.PosteriorIRF.HPDsup.' name ' = HPDIRF(:,2,j,i);']);
if max(y1) - min(y1) > 1e-10 ,
nplo=nplo+1;
if mod(nplo,nsubplo)==1,
figure('name',['Orthogonalised shocks to ',varargin{i}, ' in perc''s'])
nfig=nfig+1;
end
subplot(nn,nn,nplo)
plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
hold on
% for k = 1:9
% plot(1:options_.irf,DistribIRF(:,k,j,i),'-g','linewidth',0.5)
% end
for k = 1:2
plot(1:options_.irf,HPDIRF(:,k,j,i),'-g','linewidth',0.5)
end
plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',1)
xlim([1 options_.irf]);
hold off
title(lname{j},'interpreter','none')
end
if (mod(nplo,nsubplo)==0 | j==length(cname)) & nplo,
saveas(gcf,[DirectoryName '\' fname_,'_Bayesian_IRF_',varargin{i},'_P',int2str(nfig)])
eval(['print -depsc2 ' DirectoryName '\' fname_,'_Bayesian_IRF_',varargin{i},'_P',int2str(nfig)]);
eval(['print -dpdf ' DirectoryName,'\',fname_,'_Bayesian_IRF_',varargin{i},'_P',int2str(nfig)]);
close(gcf)
nplo=0;
end
end
end

330
matlab/gsa_sdp_fn.m
View File

@ -1,330 +0,0 @@
function gsa_ = gsa_sdp_fn(y, x, pshape, p1, p2, tvp, fname, pnames, ifig, ipred, nvr, ialias, nvrT, yi_anal)
% copyright Marco Ratto 2006
%
% function gsa_ =
% gsa_sdp_fn(y, x, pshape, p1, p2, tvp, k, pnames, ifig, ipred, nvr, ialias, nvrT, yi_anal)
%
% INPUTS
% y, the output
% x, the inputs
% pshape, input distributions
% p1, p2, parameters of input distributions
% tvp: tipe of estimation
% 0 = load previous estimation
% 1 = RW for all parameters
% 2 = IRW for all parameters
% -1 = RW for all parameters, fixed nvr's (no estimation)
% -2 = IRW for all parameters, fixed nvr's (no estimation)
% an array of length(x) specifies the process for each parameter
% fname, file name to save the analysis
% pnames: names of the params
% ifig: 1 plot figures, 0 no plots
% ipred: use polynomial interpolation to predict model output
% nvr
% yi_anal, analytical value of the first order HDMR terms (optional)
%
% OUTPUT
% gsa_.univariate.f
% gsa_.univariate.fs
% gsa_.univariate.fses
% gsa_.univariate.si
% gsa_.univariate.si_std
% gsa_.multivariate.f
% gsa_.multivariate.fs
% gsa_.multivariate.fses
% gsa_.multivariate.si
% gsa_.multivariate.si_std
% gsa_.multivariate.stat
% gsa_.x0
% gsa_.xx
% gsa_.y
%
% f, function estimates
% fs, sorted function estimates
% fses, sorted standard error of function estimates
% si, sensitivity indices
% si_std, st. error of sensitivity indices
% stat, euristic tstat for significance of f
% xx, transformed inputs used (rank-transformed)
% x0, original inputs
% y, output
[kn, np]=size(x);
if nargin<6 | isempty(tvp),
tvp=-2;
TVP=ones(1,np).*(abs(tvp)-1);
elseif length(tvp)==1,
if tvp>0
TVP=ones(1,np).*(tvp-1);
else
TVP=ones(1,np).*(abs(tvp)-1);
end
else
TVP=tvp-1;
end
if nargin<7 | isempty(fname),
fname='gsa_sdp';
end
if nargin<9 | isempty(ifig),
ifig=0;
end
if nargin<10 | isempty(ipred),
ipred=0;
end
if nargin<12 | isempty(ialias),
ialias=0;
end
if nargin<13 | isempty(nvrT),
nvrT=16;
end
if nargin<14 | isempty(yi_anal),
yi_anal=[];
end
yin=y;
[y,my,sy]=stand_(y);
x0=x;
% x=stand_(x);
if tvp~=0,
[xcum, paracdf] = cumdens(x,pshape,p1,p2);
[dum, ix]=sort(x);
[dum, iy]=sort(ix);
xx=kron([1:kn]',ones(1,np));
xx=xx(iy)./kn;
x=xx;
x=x+1;
nvr1 = ( 2*sin( 0.5*(2*pi/(kn/2)) ) )^2;
nvr2 = ( 2*sin( 0.5*(2*pi/(kn/2)) ) )^4;
%nvr2 = ( 2*sin( 0.5*(2*pi/(3*(kn/4))) ) )^4;
opts=[40 0.001 0 1 1 0]; %[iter con meth sm ALG plotopt]
for j=1:np,
lastwarn('')
if tvp<0
if tvp==-1
[fit(:,j),fitse(:,j),par(:,j),parse(:,j),zs(:,j),pars(:,j),parses(:,j),rsq(j),nvre(j)]= ...
sdr(y,x(:,j),x(:,j), 0, nvr1,[1 -1 -1 -1 -1 -1]);
else
[fit(:,j),fitse(:,j),par(:,j),parse(:,j),zs(:,j),pars(:,j),parses(:,j),rsq(j),nvre(j)]= ...
sdr(y,x(:,j),x(:,j), 1, nvr2,[1 -1 -1 -1 -1 -1]);
end
else
[fit(:,j),fitse(:,j),par(:,j),parse(:,j),zs(:,j),pars(:,j),parses(:,j),rsq(j),nvre(j)]= ...
sdr(y,x(:,j),x(:,j), TVP(j), -1,[1 -1 -1 -1 -1 -1],[],[],[],0);
end
dwarn=lastwarn;
if length(dwarn)>0,
nvre(j)=0;
end
m0=mean(zs(:,j).*pars(:,j));
vi(1, j) = var(zs(:,j).*pars(:,j));
vi(2, j) = max(abs([(par(:,j)-mean(par(:,j)))./parse(:,j)])) ;
rT(1, j) = rsq(j);
fu(:,j) =(x(:,j).*par(:,j)-m0).*sy;
fus(:,j)=zs(:,j).*pars(:,j).*sy;
fuses(:,j)=zs(:,j).*parses(:,j).*sy;
lmod(j) = {'RW'};
if TVP(j), lmod(j)={'IRW'}; end,
end
%imult = find(rT(1,:)>0.005);
[dum, imult]=sort(-vi(2,:));
%imult=[11:15];
%imult=1:np;
if isempty(imult),
mrT=max(rT(1,:));
imult = find(rT(1,:)>(0.05*mrT));
end
%[fit,fitse,par,parse,zs,pars,parses,rsq,nvre]=sdr(y,x,x,ones(1,np), ones(1,np).*(-4));
opts=[40 0.001 0 1 1 0]; %[iter con meth sm ALG plotopt]
%nvrP=ones(1,np).*(-2).*(nvre~=0);
y0=y;
if tvp<0,
if tvp==-1
nvrP=ones(1,np).*nvr1;
else
nvrP=ones(1,np).*nvr2;
end
else
nvrP=ones(1,np).*(-2);
if ialias,
nvr0=nvre.*(nvre~=0)+ones(1,np).*(1.e-10).*(nvre==0);
nvrP(ialias)=nvre(ialias);
nvr0(ialias)=nvre(ialias);
y0 = y0 - par(:,ialias).*x(:,ialias);
imult=imult(find(imult~=ialias));
else
nvr0=nvre.*(nvre~=0)+ones(1,np).*(1.e-10).*(nvre==0);
end
end
parM=par;
parseM=parse;
zsM=zs;
parsM=pars;
parsesM=parses;
nvreM=nvre;
if tvp<0
[fit2,fitse2,par2,parse2,zs2,pars2,parses2,rsq2,nvre2]= ...
sdr(y0,x(:,imult),x(:,imult),abs(tvp)-1,nvrP(imult));
else
[fit2,fitse2,par2,parse2,zs2,pars2,parses2,rsq2,nvre2]= ...
sdr(y0,x(:,imult),x(:,imult),TVP(imult),nvrP(imult),opts,[],[],nvr0(imult),0);
end
if any(nvre2>nvrT),
ik=find(nvre2>nvrT);
nvre2(ik)=ones(1,length(ik)).*nvrT;
[fit2,fitse2,par2,parse2,zs2,pars2,parses2,rsq2,nvre2]= ...
sdr(y0,x(:,imult),x(:,imult),TVP(imult),nvre2);
end
parM(:,imult)=par2;
parseM(:,imult)=parse2;
zsM(:,imult)=zs2;
parsM(:,imult)=pars2;
parsesM(:,imult)=parses2;
nvreM(imult)=nvre2;
for j=1:np,
m0=mean(zsM(:,j).*parsM(:,j));
v_di_e(1, j) = var(zsM(:,j).*parsM(:,j));
v_di_e(2,j) = mean((zsM(:,j).*parsesM(:,j)).^2);
stat(1, j) = max(abs([(parM(:,j)-mean(parM(:,j)))./parseM(:,j)])) ;
f(:,j)=(x(:,j).*parM(:,j)-m0).*sy;
fs(:,j)=(zsM(:,j).*parsM(:,j)-m0).*sy;
fses(:,j)=zsM(:,j).*parsesM(:,j).*sy;
end
pcoef = ChRegressGSA(xcum,8,1,f);
for ij=1:np,
scoef{ij}=['[',num2str(pcoef(ij,:),'%20.16e '),']'];
gg{ij}=inline(['polyval(',scoef{ij},',',char(paracdf{ij}),')'],'X');
end
gsa_.univariate.f0=my;
gsa_.univariate.f=fu;
gsa_.univariate.fs=fus;
gsa_.univariate.fses=fuses;
gsa_.univariate.nvr=nvre;
gsa_.univariate.r2=rsq2;
gsa_.univariate.si=vi(1,:);
gsa_.univariate.si_std=vi(2,:);
gsa_.multivariate.f0=my;
gsa_.multivariate.f=f;
gsa_.multivariate.fs=fs;
gsa_.multivariate.fses=fses;
gsa_.multivariate.nvr=nvreM;
gsa_.multivariate.r2=rsq2;
gsa_.multivariate.si=v_di_e(1,:);
gsa_.multivariate.si_std=v_di_e(2,:);
gsa_.multivariate.stat=stat;
gsa_.polynomial.gpred=gg;
gsa_.x0=x0;
gsa_.xx=xx;
gsa_.y=yin;
save(fname,'gsa_')
else
load(fname)
% xx=gsa_.xx;
% yin=gsa_.y;
% [xcum, paracdf] = cumdens(x,pshape,p1,p2);
% pcoef = ChRegressGSA(xcum,8,1,gsa_.multivariate.f);
% for ij=1:np,
% scoef{ij}=['[',num2str(pcoef(ij,:),'%20.16e '),']'];
% gg{ij}=inline(['polyval(',scoef{ij},',',char(paracdf{ij}),')'],'X');
% end
% gsa_.univariate.f0=my;
% gsa_.multivariate.f0=my;
% gsa_.multivariate.gpred=gg;
% save(fname,'gsa_')
end
if ifig
nplo = ceil(sqrt(np));
if nplo*(nplo-2)>=np,
nplo2=nplo-2;
elseif nplo*(nplo-1)>=np,
nplo2=nplo-1;
elseif (nplo+1)*(nplo-1)>=np,
nplo2=nplo-1;
nplo=nplo+1;
else
nplo2=nplo;
end
if nargin<8 | isempty(pnames),
pnames=['X_{',num2str(1),'}'];
for j=2:np,
pnames=str2mat(pnames,['X_{',num2str(j),'}']);
end
end
figure('name',['Nrun=',num2str(kn),', ',fname,', SDP multivariate']),
f=gsa_.multivariate.f;
fs=gsa_.multivariate.fs;
fses=gsa_.multivariate.fses;
stat=gsa_.multivariate.stat;
v_di_e(1,:)=gsa_.multivariate.si;
for j=1:np,
%subplot(3,4,j)
subplot(nplo2,nplo,j)
xp=sort(gsa_.x0(:,j));
%xp=zsM(:,jj)-1;
if exist('yi_anal') & ~isempty(yi_anal),
%plot(xp, yi_anal(:,j)./sy,'r:'), hold on,
plot(xp, yi_anal(:,j),'r:'), hold on,
end
plot(xp, fs(:,j),'k', xp, 3.09.*fses(:,j),':k', xp, -3.09.*fses(:,j),':k')
set(gca,'xlim',[min(xp) max(xp)])
lmod(j) = {'I1'};
if TVP(j), lmod(j)={'I2'}; end,
title([deblank(pnames(j,:)),', Si=',num2str(v_di_e(1,j),'%5.2f')],'fontsize',9,'interpreter','none')
end
saveas(gcf,[fname,'_SE_multiv'])
eval(['print -depsc2 ' fname '_SE_multiv']);
eval(['print -dpdf ' fname '_SA_multiv']);
close(gcf)
% figure('name',['Nrun=',num2str(kn),', ',fname,', SDP univariate']),
% f=gsa_.univariate.f;
% fs=gsa_.univariate.fs;
% fses=gsa_.univariate.fses;
% vi(1,:)=gsa_.univariate.si;
% vi(2,:)=gsa_.univariate.si_std;
% for j=1:np,
% subplot(nplo2,nplo,j)
% xp=sort(x0(:,j));
% if exist('yi_anal') & ~isempty(yi_anal),
% plot(xp, yi_anal(:,j),'r'), hold on,
% end
% plot(xp, [fs(:,j) 3.09.*fses(:,j) -3.09.*fses(:,j)])
% set(gca,'xlim',[min(xp) max(xp)])
% lmod(j) = {'I1'};
% if TVP(j), lmod(j)={'I2'}; end,
% title([deblank(pnames(j,:)),', ',lmod{j},', tstat=',num2str(vi(2,j),'%5.2f'), ...
% ', S=',num2str(vi(1,j),'%5.2f')])
% end
% %saveas(gcf,[fname,'_SE_univ'])
% close(gcf)
end
%end
if ipred,
gg=gsa_.polynomial.gpred;
for j=1:np;
ff(:,j)=gg{j}(x0(:,j));
end
gsa_.polynomial.ypred=sum(ff,2)+gsa_.multivariate.f0;
gsa_.polynomial.r2=1-cov(gsa_.polynomial.ypred-yin)/cov(yin);
for j=1:12; subplot(3,4,j), plot(x0(:,j),ff(:,j),'.',gsa_.x0(:,j),f(:,j),'b.'), end
gsa_.xpred=x0;
% nvreM=gsa_.multivariate.nvr;
% x=[gsa_.x0; x0];
% xcum = cumdens(x,pshape,p1,p2);
% [yin,my,sy]=stand_(gsa_.y);
% yf=[yin; NaN.*ones(size(x0,1),1)];
% yfit2= sdp(yf,xcum,xcum,TVP,nvreM);
% gsa_.ypred_sdr = yfit2(length(yin)+1:end).*sy+my;
save(fname,'gsa_')
end

235
matlab/ksirw2c.m
View File

@ -1,235 +0,0 @@
function [Xkk,Pkk,er,ykk1,Xkk1,yp,ep,ykk,XkN,ys,PkN,ers]=...
ksirw2c(X,RWtype,Miss,Interv,F,P0,Q,iout,x0,nf,intD,QPflag)
% [Xkk,Pkk,er,ykk1,Xkk1,yp,ep,ykk,XkN,ys,PkN,ers]=...
% ksirw2c(X,RWtype,Miss,Interv,F,P0,Q,iout,x0,nf,intD,QPflag)
% 1 2 3 4 5 6 7 8 9 10 11 12
% Fast version of KALMSMO for use with SDR
% No interventions, no missing values, only RW and IRW single parameter models
% Execution time gain - up to 10 times.
% WT, September 2004
% Input:
% X: (Ndat,Np+1) - data (y-s in the first column)
% RWtype (Np) - number of aux. parameters (eg. slopes)
% for given parameters
% Miss: (Ndat,1) - missing data flags (1 - missing sample)
% Interv: (Ndat,1) - intervention flags (1 - intervention)
% F : (sqare matrix Np2) - transition matrix
% Q : (sqare matrix Np2) - noise variance (nvr) matrix
% P0 : (sqare matrix Np2) - parameters cov. matrix (Pkk) - initial
% iout (integer) - intermediate results on/off
% x0 - initial cond. for KF
% nf - max forecasting horizon
% intD - diagonal of Q matrix defining the
% intervention type (size of diag(Q))
% Output:
% XkN: as Xkk but smoothed
% ys : smoothed output
% PkN: (square matrix Np2 repeated 1:Ndat) - smoothed par. cov. mtx
% Xkk: (Ndat,Np2) - filtered parameters history
% Pkk : - P(k|k) history
% er : (Ndat) - normalised error
% ykk1: (Ndat) - one step ahead prediction
% yp: (Ndat,nf-1) - predictions: yp(k,j)=y(k+j|k-1)
% ep: (Ndat,nf-1) - nrm. fact.of prediction errors:
% - ep(k|j)=s(k+j|k-1)
% ers : (Ndat) - normalised error (smoothing)
% Q and P algorithms implemented (QPflag==1 : Q alg)
% Copyright (c) 2004 by CRES, Lancaster University, United Kingdom
% Authors : Peter Young, Wlodek Tych, Diego Pedregal, James Taylor
% Revision history :
% 20/09/2004, WT, converted to ksrivc from kalmsmo
% simplifications:
% RWtype = 1 or 0
% F = [1 1;0 1] or 1
% Q = [0 0;0 nvr] or nvr
% no interventions
% no missing data
% no forecasting built in
yp=[];ep=[];ykk=[];XkN=[];ys=[];PkN=[];ers=[];
% sizes and initial conds.
[Ndat Np]=size(X);Np=Np-1;Np2=length(F);
P=P0;
% defaults
if nargin<12
QPflag=1;
end
if (nargin<9); x0=zeros(Np2,1); end
if isempty(x0),x0=zeros(Np2,1); end
if nargin<8, iout=0; end % iout
if isempty(iout), iout=1; end
er=zeros(Ndat,1);
perr=zeros(Ndat,1);
Xkk=zeros(Np2,Ndat); Xkk1=[]; % Xkk1=Xkk;
Pkk=zeros(Np2,Np2*Ndat); Pkk(:,1:Np2)=P0;
ykk1=zeros(Ndat,1); ykk=ykk1;
PHmat=zeros(Np2,Ndat);
nvr=Q(Np2,Np2);
Xhat=zeros(Np2,1);
if Np2==2 % **************************** IRW
if ~isempty(x0), Xhat(1)=x0(1);Xhat(2)=x0(2);end
for it=1:Ndat
u=X(it,2);
% KF prediction step
Xhat=[ Xhat(1)+Xhat(2); Xhat(2)]; % state X(k|k-1)
P=[ P(1,1)+P(2,1)+P(1,2)+P(2,2) P(1,2)+P(2,2);P(2,1)+P(2,2) P(2,2)+nvr]; % cov. matrix
HP=[ u*P(1,1), u*P(1,2)];
div=1+u*P(1,1)*u;
er(it)=div;div=1/div;
yhat=u*Xhat(1);ykk1(it)=yhat; % calc. and save save y(k|k-1)
perr(it)=X(it,1)-yhat; % innovation
% KF correction step
Xhat=Xhat+HP'*perr(it)*div; % X(k|k)
yhat=u*Xhat(1); % yhat(k|k)
P=P-u*u*[P(1,1)*P(1,1) P(1,1)*P(1,2);P(1,2)*P(1,1) P(1,2)*P(1,2)]*div;
ykk(it)=yhat; % save
Xkk(:,it)=Xhat;
Pkk(:,(it-1)*Np2+1:it*Np2)=P;
PHmat(:,it)=u*[ P(1,1);P(2,1)];
end; % of filtering loop
L=zeros(2,1);
XkN=zeros(Np2,Ndat); XkN(:,Ndat)=Xhat;
ys=zeros(Ndat,1); ys(Ndat)=Xhat(1);
PkN=zeros(Np2,Np2*Ndat);
PkN(:,(Ndat-1)*Np2+1:Ndat*Np2)=P;
ers=zeros(Ndat,1);
sP=P;
ers(Ndat)=er(Ndat);
% smoothing loop
for it=Ndat:-1:1
u=X(it,2);
XkN(:,it)=Xhat;
if QPflag,
ys(it)=u*Xhat(1);
else
pkk=Pkk(:,(it-1)*Np2+1:(it)*Np2);
Xhat=Xkk(:,it)-[(pkk(1,1)+pkk(1,2))*L(1)+pkk(1,2)*L(2); (pkk(2,1)+pkk(2,2))*L(1)+pkk(2,2)*L(2)];
end
% if Miss(it)==0
L=[(1-PHmat(1,it)*u)*(L(1)-perr(it)*u)-PHmat(2,it)*u*(L(1)+L(2)); L(1)+L(2)];
% else
% L=[L(1);L(1)+L(2)];
% end
if QPflag,
Xhat=[Xhat(1)-Xhat(2)-nvr*L(2);Xhat(2)+nvr*L(2)];
else
ys(it)=u*Xhat(1);
XkN(:,it)=Xhat;
end
if (it<Ndat)&(nargout>10) % calculation of PkN requested
P=Pkk(:,(it-1)*Np2+1:it*Np2);
Pdd=(P(1,1)*P(2,2)+P(1,1)*nvr+P(2,1)*nvr+P(1,2)*nvr+P(2,2)*nvr-P(2,1)*P(1,2));
Pcc=(-P(1,1)*P(2,2)+P(2,1)*P(1,2));
Paa=(Pcc/Pdd);
Pbb=(sP(1,1)-P(1,1)-P(2,1)-P(1,2)-P(2,2));
Pff=(sP(2,1)-P(2,1)-P(2,2));
Pgg=(nvr*(P(2,1)+P(2,2))/Pdd*(sP(1,2)-P(1,2)-P(2,2))-Paa*(sP(2,2)-P(2,2)-nvr));
Phh=(P(1,1)*P(2,2)+P(1,1)*nvr+P(1,2)*nvr-P(2,1)*P(1,2));
Pbd=Pbb/Pdd;
Phd=(Phh/Pdd);
Pah=(sP(2,2)-P(2,2)-nvr);
Pah1=Phd*(sP(1,2)-P(1,2)-P(2,2));
Ph0=(Phh*Pbd+Paa*Pff);
Pn=nvr*(P(2,1)+P(2,2));
sP=[P(1,1)+Ph0*Phd+(Pah1+Paa*Pah)*Paa P(1,2)+Ph0*(Pn/Pdd)-(Pah1+Paa*Pah)*Paa;
P(2,1)+(Pn*Pbd-Paa*Pff)*Phd+Pgg*Paa P(2,2)+(Pn*Pbd-Paa*Pff)*(Pn/Pdd)-Pgg*Paa];
PkN(:,(it-1)*Np2+1:it*Np2)=sP;
ers(it)=1+u*u*sP(1,1);
end
if(abs(iout)>2)&~rem(it,10) % display current est.
disp(Xhat(ipp)'); sPd=diag(sP)';
disp(sqrt(sPd(ipp)))
end
end;
ys=ys'; % consistent with previous versions
else % if Np2==2 ************************************ RW
if ~isempty(x0), Xhat=x0;end
for it=1:Ndat
u=X(it,2);
Xhat=Xhat;
P=P+nvr;
HP=u*P;
div=1+u*P*u;
er(it)=div;div=1/div;
yhat=u*Xhat;ykk1(it)=yhat; % calc. and save save y(k|k-1)
perr(it)=X(it,1)-yhat;
Xhat=Xhat+HP*perr(it)*div;
yhat=u*Xhat;
P=P-u*u*P*P*div;
ykk(it)=yhat; % save
Xkk(:,it)=Xhat;
Pkk(it)=P;
PHmat(it)=u*P;
end; % of filtering loop
L=0;
XkN=zeros(Np2,Ndat); XkN(:,Ndat)=Xhat;
ys=zeros(Ndat,1); ys(Ndat)=Xhat(1);
PkN=zeros(Np2,Np2*Ndat);
PkN(:,(Ndat-1)*Np2+1:Ndat*Np2)=P;
ers=zeros(Ndat,1);
sP=P;
ers(Ndat)=er(Ndat);
for it=Ndat:-1:1
u=X(it,2);
XkN(:,it)=Xhat;
if QPflag,
ys(it)=u*Xhat;
else
pkk=Pkk(:,(it-1)*Np2+1:(it)*Np2);
Xhat=Xkk(it)-Pkk(it)*L;
end
% if Miss(it)==0
L=(1-u*PHmat(it))*(L-u*perr(it));
% else
% L=L;
% end
if QPflag,
Xhat=Xhat+nvr*L;
else
ys(it)=u*Xhat;
XkN(it)=Xhat;
end
if (it<Ndat)&(nargout>10)
P=Pkk(it);
Pp=P+Q;
PT=P/Pp;
sP=P+PT*(sP-Pp)*PT;
PkN(it)=sP;
ers(it)=1+sP;
end
end;
ys=ys'; % consistent with previous versions
end % if Np2==2
% end of m-file

BIN
matlab/ksirw2c.p
View File

Binary file not shown.

19
matlab/mh_moma.m
View File

@ -1,19 +0,0 @@
function th_std = mh_moma(var_list_)
global M_ oo_
DirectoryName = CheckPath('Metropolis');
load([DirectoryName,'\',M_.fname,'_param_irf1.mat'])
h = waitbar(0,'Computing MH theoretical moments...');
B=size(stock,1);
th_std=zeros(size(var_list_,1),B);
for j=1:B,
set_all_parameters(stock(j,:))
stoch_simul(var_list_);
th_std(:,j)=sqrt(diag(oo_.var));
waitbar(j/B,h)
end
close(h)

30
matlab/myDefault.m
View File

@ -1,30 +0,0 @@
function myDefault
set(0,'DefaultFigurePaperType','a4')
%set(0,'DefaultFigurePaperOrientation','landscape')
set(0,'DefaultFigureColor',[0 0 0])
set(0,'DefaultTextColor',[1 1 1])
set(0,'DefaultAxesUnits','normalized')
set(0,'DefaultTextUnits','normalized')
set(0,'DefaultAxesLineStyleOrder', ...
{'-',':','--','-.'})
set(0,'DefaultAxesColor',[0 0 0],...
'DefaultAxesXColor', [1 1 1], ...
'DefaultAxesYColor', [1 1 1], ...
'DefaultAxesZColor', [1 1 1])
set(0, ...
'DefaultAxesColorOrder', ...
[1 1 0; ... %yellow
1 1 1; ... % white
0 1 1; ... % cyan
1 0 1; ... % magenta
1 0 0; ... % red
0 1 0; ... % green
0 0 1]) % blue
%set(0,'DefaultPatchmarkeredgecolor',[1 1 1])
set(0,'DefaultPatchEdgeColor',[1 1 1])
%map = colormap;
%colormap(map([end:-1:1],:));
%colormap(hsv);

8
matlab/name2cell.m
View File

@ -1,8 +0,0 @@
function name2cell(varargin);
cname = varargin(2:end);
if strcmp(varargin{1}(1),'-'),
assignin('caller',varargin{1}(2:end),cname)
else
disp('The name of the cell variable MUST start with - (minus)!!!')
end

61
matlab/plot_IRF.m
View File

@ -1,61 +0,0 @@
function plot_IRF(xname, vargin, s1, texname)
global M_ oo_ options_
fname_ =M_.fname;
DirectoryName = CheckPath('Output');
if nargin<4,
for j=1:M_.endo_nbr,
texname{j}=deblank(M_.endo_names_tex(j,:));
end
end
iendo=[];
for i=1:length(vargin)
iendo=[iendo, strmatch(vargin{i},M_.endo_names,'exact')];
end
iexo=[];
for i=1:length(xname)
iexo=[iexo, strmatch(xname{i},M_.exo_names,'exact')];
end
for i=1:length(iexo),
nfig=0;
nplo=0;
for j=1:length(vargin),
name = [vargin{j} '_' deblank(M_.exo_names(iexo(i),:))];
eval(['MeanIRF=s1.' name,';']);
if max(abs(MeanIRF)) > 1e-6 ,
nplo=nplo+1;
if mod(nplo,9)==1,
figure('name',['Orthogonalised shocks to ',deblank(M_.exo_names(iexo(i),:))])
nfig=nfig+1;
end
subplot(3,3,nplo)
plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
hold on,
plot(1:options_.irf,MeanIRF(:),'-k','linewidth',1)
xlim([1 options_.irf]);
hold off
if options_.TeX,
if nargin<4,
title(texname{iendo(j)},'interpreter','tex')
else
title(texname{j},'interpreter','tex')
end
else
title(vargin{j},'interpreter','none')
end
end
if (mod(nplo,9)==0 | j==length(vargin)) & nplo,
saveas(gcf,[DirectoryName '\' fname_,'_IRF_',deblank(M_.exo_names(iexo(i),:)),'_',int2str(nfig)])
eval(['print -depsc2 ' DirectoryName '\' fname_,'_IRF_',deblank(M_.exo_names(iexo(i),:)),'_',int2str(nfig)]);
eval(['print -dpdf ' DirectoryName,'\',fname_,'_IRF_',deblank(M_.exo_names(iexo(i),:)),'_',int2str(nfig)]);
close(gcf)
nplo=0;
end
end
end

63
matlab/plot_IRF_post.m
View File

@ -1,63 +0,0 @@
function plot_IRF_post(xname, vargin, texname)
global M_ oo_ options_
fname_ =M_.fname;
DirectoryName = CheckPath('Output');
iendo=[];
for i=1:length(vargin)
iendo=[iendo, strmatch(vargin{i},M_.endo_names,'exact')];
end
if nargin<3,
for j=1:length(iendo),
texname{j}=deblank(M_.endo_names_tex(iendo(j),:));
end
end
iexo=[];
for i=1:length(xname)
iexo=[iexo, strmatch(xname{i},M_.exo_names,'exact')];
end
for i=1:length(iexo),
nfig=0;
nplo=0;
for j=1:length(vargin),
name = [vargin{j} '_' deblank(M_.exo_names(iexo(i),:))];
eval(['MeanIRF=oo_.PosteriorIRF.Mean.' name,';']);
eval(['DistribIRF=oo_.PosteriorIRF.Distribution.' name ';']);
eval(['HPDsup=oo_.PosteriorIRF.HPDsup.' name ';']);
eval(['HPDinf=oo_.PosteriorIRF.HPDinf.' name ';']);
if max(abs(MeanIRF)) > 1e-6 ,
nplo=nplo+1;
if mod(nplo,9)==1,
figure('name',['Orthogonalised shocks to ',deblank(M_.exo_names(iexo(i),:))])
nfig=nfig+1;
end
subplot(3,3,nplo)
patch([[1:options_.irf] [options_.irf:-1:1]],[HPDsup' HPDinf(end:-1:1)'],[0.75 0.75 0.75])
hold on
plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
% for k = 1:9
% plot(1:options_.irf,DistribIRF(:,k),'-g','linewidth',0.5)
% end
plot(1:options_.irf,MeanIRF(:),'-k','linewidth',1)
xlim([1 options_.irf]);
hold off
if options_.TeX,
title(texname{j},'interpreter','tex')
else
title(vargin{j},'interpreter','none')
end
end
if (mod(nplo,9)==0 | j==length(vargin)) & nplo,
saveas(gcf,[DirectoryName '\' fname_,'_Bayesian_IRF_',deblank(M_.exo_names(iexo(i),:)),'_',int2str(nfig)])
eval(['print -depsc2 ' DirectoryName '\' fname_,'_Bayesian_IRF_',deblank(M_.exo_names(iexo(i),:)),'_',int2str(nfig)]);
eval(['print -dpdf ' DirectoryName,'\',fname_,'_Bayesian_IRF_',deblank(M_.exo_names(iexo(i),:)),'_',int2str(nfig)]);
close(gcf)
nplo=0;
end
end
end

37
matlab/plot_chain.m
View File

@ -1,37 +0,0 @@
function plot_chain
global M_ options_ estim_params_ bayestopt_
npar = estim_params_.np+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.nvx;
nblck = options_.mh_nblck;
DirectoryName = CheckPath('metropolis');
OutDirectoryName = CheckPath('output');
load([DirectoryName '/' M_.fname '_mh_history'])
TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
ifig=0;
for j=1:npar,
if mod(j,9)==1,
ifig=ifig+1;
iplot=0;
hh=figure('Name',['Markov Chain plot ',int2str(ifig)]);
end
iplot=iplot+1;
Draws = GetAllPosteriorDraws(j,1,1,TotalNumberOfMhFiles,TotalNumberOfMhDraws);
Draws=reshape(Draws,TotalNumberOfMhDraws,nblck);
Draws = stand_(Draws);
subplot(3,3,iplot)
plot(Draws)
hold on, plot(cumsum(Draws)./[1:length(Draws)]','r')
title(bayestopt_.name{j},'interpreter','none')
if mod(j,9)==0 | j==npar,
saveas(hh,[OutDirectoryName '\' M_.fname '_chain_' int2str(ifig)])
eval(['print -depsc2 ' OutDirectoryName '/' M_.fname '_chain_' int2str(ifig)]);
eval(['print -dpdf ' OutDirectoryName '/' M_.fname '_chain_' int2str(ifig)]);
close(hh)
end
end

78
matlab/plot_conditionalIRF_post.m
View File

@ -1,78 +0,0 @@
function plot_conditionalIRF_post(xname, vargin, texname)
global M_ oo_ options_ bayestopt_
fname_ =M_.fname;
DirectoryName = CheckPath('Output');
iendo=[];
for i=1:length(vargin)
iendo=[iendo, strmatch(vargin{i},M_.endo_names,'exact')];
end
if nargin<3,
for j=1:length(iendo),
texname{j}=deblank(M_.endo_names_tex(iendo(j),:));
end
end
iexo=[];
for i=1:length(xname)
iexo=[iexo, strmatch(xname{i},M_.exo_names,'exact')];
end
y=getIRFRuns(iendo,iexo);
DirName = CheckPath('Metropolis');
load([DirName,'\',M_.fname,'_param_irf1.mat'])
for i=1:length(iexo),
nfig=0;
nplo=0;
ipar=strmatch(xname{i},bayestopt_.name,'exact');
for j=1:length(vargin),
y0=squeeze(y(:,j,i,:));
if ~isempty(ipar),
for k=1:size(y0,1),
x0=stand(stock(:,ipar));
b=regress(y0(k,:)',[ones(500,1) x0]);
y0(k,:)=y0(k,:)-x0'.*b(2);
end
end
for k=1:options_.irf,
[MeanIRF(k),MedianIRF(k),VarIRF(k),HPDIRF(k,:),DistribIRF(k,:)] = ...
posterior_moments(y0(k,:),0);
end
name = [vargin{j} '_' deblank(M_.exo_names(iexo(i),:))];
if max(abs(MeanIRF)) > 1e-6 ,
nplo=nplo+1;
if mod(nplo,9)==1,
figure('name',['Orthogonalised shocks to ',deblank(M_.exo_names(iexo(i),:))])
nfig=nfig+1;
end
HPDinf = HPDIRF(:,1);
HPDsup = HPDIRF(:,2);
subplot(3,3,nplo)
patch([[1:options_.irf] [options_.irf:-1:1]],[HPDsup' HPDinf(end:-1:1)'],[0.75 0.75 0.75])
hold on
plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
% for k = 1:9
% plot(1:options_.irf,DistribIRF(:,k),'-g','linewidth',0.5)
% end
plot(1:options_.irf,MeanIRF,'-k','linewidth',1)
xlim([1 options_.irf]);
hold off
if options_.TeX,
title(texname{j},'interpreter','tex')
else
title(vargin{j},'interpreter','none')
end
end
if (mod(nplo,9)==0 | j==length(vargin)) & nplo,
saveas(gcf,[DirectoryName '\' fname_,'_Bayesian_cond_IRF_',deblank(M_.exo_names(iexo(i),:)),'_',int2str(nfig)])
eval(['print -depsc2 ' DirectoryName '\' fname_,'_Bayesian_cond_IRF_',deblank(M_.exo_names(iexo(i),:)),'_',int2str(nfig)]);
eval(['print -dpdf ' DirectoryName,'\',fname_,'_Bayesian_cond_IRF_',deblank(M_.exo_names(iexo(i),:)),'_',int2str(nfig)]);
close(gcf)
nplo=0;
end
end
end

81
matlab/plot_fcast.m
View File

@ -1,81 +0,0 @@
function [fcast] = plot_fcast(k_, varargin);
global options_ bayestopt_ M_ oo_
ys_ = oo_.steady_state;
lgy_ = M_.endo_names;
fname_ = M_.fname;
if isempty(k_), k_=options_.nk; end
if ischar(k_),
k_=eval(k_);
end
if exist(options_.datafile)
instr = options_.datafile;
else
instr = ['load ' options_.datafile];
end
eval(instr);
if ~exist('T'),
temp = eval(deblank(options_.varobs(1,:)));
T=[1:length(temp)]';
clear temp;
end
if isempty(varargin)
lgobs_= options_.varobs;
mfys=bayestopt_.mfys;
else
mfys=[];
for j=1:length(varargin),
dum = strmatch(varargin{j},lgy_,'exact');
mfys = [mfys dum];
if j==1,
lgobs_ = varargin{j};
else
lgobs_ = str2mat(lgobs_,varargin{j});
end
end
end
nobs = size(lgobs_,1);
fobs = options_.first_obs;
if options_.loglinear == 1
constant = log(ys_(mfys));
else
constant = ys_(mfys);
end
lst_obs = options_.nobs+options_.first_obs - 1;
%disp(' ')
%disp(' ')
%disp(' k-step RMSE')
for ifig = 1:ceil(size(lgobs_,1)/9),
figure('Name',['DSGE ',int2str(k_),'-step forecasts']),
for j=1+9*(ifig-1):min(9*ifig,size(lgobs_,1)),
subplot(3,3,j-9*(ifig-1)),
vj=deblank(lgobs_(j,:));
i = strmatch(vj,lgy_,'exact');
for l=1:k_
eval(['vp(l)= oo_.KStepPredictions.', ...
vj,'(options_.nobs+l,l)+constant(j);'])
end
plot([T(lst_obs)+0.25:0.25:T(lst_obs)+k_/4], vp)
hold on,
if exist(vj,'var')
eval(['vd= ',vj,';'])
plot([T(lst_obs-8:end)], vd(lst_obs-8:end),'r')
else
eval(['vd= oo_.SmoothedVariables.',vj,'+constant(j);'])
plot(T(lst_obs-8:end), vd(end-8:end) ,'r')
end
title(vj,'interpreter','none')
hh=get(gca,'children');
end
saveas(gcf,[fname_,'_',int2str(k_),'_fcast',int2str(ifig)])
if options_.nograph
close(gcf)
end
end

100
matlab/plot_fit.m
View File

@ -1,100 +0,0 @@
function [rmse, lgobs_, vv] = plot_fit(varargin);
%function [rmse, lgobs_] = plot_fit(varargin);
global options_ bayestopt_ M_ oo_
ys_ = oo_.steady_state;
lgy_ = M_.endo_names;
fname_ = M_.fname;
texname = M_.endo_names_tex;
if exist(options_.datafile)
instr = options_.datafile;
else
instr = ['load ' options_.datafile];
end
eval(instr);
if ~exist('T','var'),
temp = eval(deblank(options_.varobs(1,:)));
T=[1:length(temp)]';
clear temp;
end
istart = options_.presample+1;
if isempty(varargin),
lgobs_= options_.varobs;
mfys=bayestopt_.mfys;
else
mfys=[];
for j=1:length(varargin),
dum = strmatch(varargin{j},lgy_,'exact');
mfys = [mfys dum];
if j==1,
lgobs_ = varargin{j};
else
lgobs_ = str2mat(lgobs_,varargin{j});
end
end
end
nobs = size(lgobs_,1);
fobs = options_.first_obs;
if options_.loglinear == 1
constant = log(ys_(mfys));
else
constant = ys_(mfys);
end
trend = constant*ones(1,options_.nobs);
%disp(' ')
%disp(' ')
%disp(' RMSE')
for ifig = 1:ceil(size(lgobs_,1)/9),
h1 = figure('Name',['DSGE 1-step ahead prediction']);
h2 = figure('Name',['DSGE Innovations']);
for j=1+9*(ifig-1):min(9*ifig,size(lgobs_,1)),
figure(h1)
subplot(3,3,j-9*(ifig-1)),
vj = deblank(lgobs_(j,:));
i = strmatch(vj,lgy_,'exact');
if ~isempty(texname),
if strmatch('\log',texname(i,:));
texname(i,1:end-1)=texname(i,2:end);
end
end
eval(['plot(T(fobs+istart-1:fobs+options_.nobs-1), [oo_.FilteredVariables.', ...
vj,'(istart:end-1)+trend(j,istart:end)'' oo_.SmoothedVariables.',...
vj,'(istart:end)+trend(j,istart:end)''])'])
if options_.TeX,
title(deblank(texname(i,:)),'interpreter','tex')
else
title(vj,'interpreter','none')
end
eval(['rmse(j) = sqrt(mean((oo_.SmoothedVariables.',vj,'(istart:end)-oo_.FilteredVariables.',vj,'(istart:end-1)).^2));'])
eval(['vv(:,j) = (oo_.SmoothedVariables.',vj,'(istart:end)-oo_.FilteredVariables.',vj,'(istart:end-1));'])
hh=get(gca,'children');
set(hh(1),'linestyle',':','color',[0 0 0])
set(hh(2),'linestyle','-','color',[0 0 0])
set(gca,'xlim',[T(fobs)-1 T(fobs+options_.nobs-1)+1])
%disp([lgobs_(j,:), sprintf('%15.5g',[rmse(j)'])])
figure(h2)
subplot(3,3,j-9*(ifig-1)),
plot(T(fobs+istart-1:fobs+options_.nobs-1), vv(:,j),'k')
if options_.TeX,
title(deblank(texname(i,:)),'interpreter','tex')
else
title(vj,'interpreter','none')
end
set(gca,'xlim',[T(fobs)-1 T(fobs+options_.nobs-1)+1])
end
saveas(h1,[fname_,'_Fit',int2str(ifig)])
eval(['print -depsc2 ' fname_,'_Fit',int2str(ifig)]);
eval(['print -dpdf ' fname_,'_Fit',int2str(ifig)]);
saveas(h2,[fname_,'_Innovations',int2str(ifig)])
eval(['print -depsc2 ' fname_,'_Innovations',int2str(ifig)]);
eval(['print -dpdf ' fname_,'_Innovations',int2str(ifig)]);
if options_.nograph
close(h1)
close(h2)
end
end

107
matlab/plot_fit_post.m
View File

@ -1,107 +0,0 @@
function plot_fit_post(varargin);
%function [rmse, lgobs_] = plot_fit(varargin);
global options_ bayestopt_ M_ oo_
dirname = CheckPath('Output');
ys_ = oo_.steady_state;
lgy_ = M_.endo_names;
fname_ = M_.fname;
if exist(options_.datafile)
instr = options_.datafile;
else
instr = ['load ' options_.datafile];
end
eval(instr);
if ~exist('T','var'),
temp = eval(deblank(options_.varobs(1,:)));
T=[1:length(temp)]';
clear temp;
end
istart = options_.presample+1;
if isempty(varargin),
lgobs_= options_.varobs;
mfys=bayestopt_.mfys;
else
mfys=[];
for j=1:length(varargin),
dum = strmatch(varargin{j},lgy_,'exact');
mfys = [mfys dum];
if j==1,
lgobs_ = varargin{j};
else
lgobs_ = str2mat(lgobs_,varargin{j});
end
end
end
nobs = size(lgobs_,1);
fobs = options_.first_obs;
if options_.loglinear == 1
constant = log(ys_(mfys));
else
constant = ys_(mfys);
end
trend = constant*ones(1,options_.nobs);
%disp(' ')
%disp(' ')
%disp(' RMSE')
y00=getFilterRuns(lgobs_);
if isempty(y00), return, end,
x = T(fobs+istart-1:fobs+options_.nobs-1);
for ifig = 1:ceil(size(lgobs_,1)/9),
h1 = figure('Name',['Bayesian DSGE 1-step ahead prediction']);
h2 = figure('Name',['Bayesian DSGE Innovations']);
for j=1+9*(ifig-1):min(9*ifig,size(lgobs_,1)),
figure(h1)
subplot(3,3,j-9*(ifig-1)),
vj = deblank(lgobs_(j,:));
i = strmatch(vj,lgy_,'exact');
for k=1:options_.nobs,
[y0(k,1),Median(k,1),Var(k,1),HPD(k,:),Distrib(k,:)] = ...
posterior_moments(y00(j,k,:),0);
end
yobs=eval(vj);
% patch([x, x(end:-1:1)], [HPD(istart:end ,1)', HPD(end:-1:istart,2)'],[0.75 0.75 0.75]);
% hold on,
plot(x, [y0(istart:end) yobs(fobs+istart-1:fobs+options_.nobs-1)])
if options_.TeX,
tnam=deblank(M_.endo_names_tex(i,:));
if strmatch('\log',tnam),
tnam=tnam(2:end);
end
title(tnam,'interpreter','tex')
else
title(vj,'interpreter','none')
end
vv(:,j) = y0(istart:end)-yobs(fobs+istart-1:fobs+options_.nobs-1);
hh=get(gca,'children');
set(hh(1),'linestyle',':','color',[0 0 0])
set(hh(2),'linestyle','-','color',[0 0 0])
set(gca,'xlim',[T(fobs)-1 T(fobs+options_.nobs-1)+1])
hold off,
%disp([lgobs_(j,:), sprintf('%15.5g',[rmse(j)'])])
figure(h2)
subplot(3,3,j-9*(ifig-1)),
plot(x, vv(:,j),'k')
if options_.TeX,
title(tnam,'interpreter','tex')
else
title(vj,'interpreter','none')
end
set(gca,'xlim',[T(fobs)-1 T(fobs+options_.nobs-1)+1])
end
saveas(h1,[dirname,'\',fname_,'_Bayesian_Fit',int2str(ifig)])
eval(['print -depsc2 ' dirname,'\',fname_,'_Bayesian_Fit',int2str(ifig)]);
eval(['print -dpdf ' dirname,'\',fname_,'_Bayesian_Fit',int2str(ifig)]);
saveas(h2,[dirname,'\',fname_,'_Bayesian_Innovations',int2str(ifig)])
if options_.nograph
close(h1)
close(h2)
end
end

68
matlab/plot_kstep.m
View File

@ -1,68 +0,0 @@
function [rmseK, lgobs_] = plot_kstep(k_, varargin);
global options_ bayestopt_ M_ oo_
ys_ = oo_.steady_state;
lgy_ = M_.endo_names;
fname_ = M_.fname;
if isempty(k_), k_=options_.nk; end
if ischar(k_),
k_=eval(k_);
end
if exist(options_.datafile)
instr = options_.datafile;
else
instr = ['load ' options_.datafile];
end
eval(instr);
if ~exist('T'),
temp = eval(deblank(options_.varobs(1,:)));
T=[1:length(temp)]';
clear temp;
end
lgobs_= options_.varobs;
mfys=bayestopt_.mfys;
for j=1:length(varargin),
dum = strmatch(varargin{j},lgy_,'exact');
mfys = [mfys dum];
lgobs_ = str2mat(lgobs_,varargin{j});
end
nobs = size(lgobs_,1);
fobs = options_.first_obs;
if options_.loglinear == 1
constant = log(ys_(mfys));
else
constant = ys_(mfys);
end
trend = constant*ones(1,options_.nobs);
istart=k_+options_.presample;
%disp(' ')
%disp(' ')
%disp(' k-step RMSE')
for ifig = 1:ceil(size(lgobs_,1)/9),
figure('Name',['DSGE ',int2str(k_),'-step ahead prediction']),
for j=1+9*(ifig-1):min(9*ifig,length(lgobs_)),
subplot(3,3,j-9*(ifig-1)),
vj=deblank(lgobs_(j,:));
i = strmatch(vj,lgy_,'exact');
eval(['plot(T(fobs+istart-1:fobs+options_.nobs-1), [oo_.KStepPredictions.', ...
vj,'(istart:end-options_.nk,k_)+trend(j,istart:end)'' oo_.SmoothedVariables.',...
vj,'(istart:end)+trend(j,istart:end)''])'])
title(vj,'interpreter','none')
hh=get(gca,'children');
set(hh(1),'linestyle',':','color',[0 0 0])
set(hh(2),'linestyle','-','color',[0 0 0])
set(gca,'xlim',[T(fobs)-1 T(fobs+options_.nobs-1)+1])
eval(['rmseK(j) = sqrt(mean((oo_.SmoothedVariables.',vj,'(istart:end)-oo_.KStepPredictions.',vj,'(istart:end-options_.nk,k_)).^2));'])
%disp([lgobs_(j,:), sprintf('%15.5g',[rmseK(j)'])])
end
saveas(gcf,[fname_,'_',int2str(k_),'_pred',int2str(ifig)])
if options_.nograph
close(gcf)
end
end

65
matlab/plot_smooth.m
View File

@ -1,65 +0,0 @@
function plot_smooth(varargin)
global options_ M_ oo_
ys_ = oo_.steady_state;
lgy_ = M_.endo_names;
fname_ = M_.fname;
texname=M_.endo_names_tex;
ifig0=0;
if strmatch('-',varargin{end})
optn = varargin{end};
varargin=varargin(1:end-1);
if strmatch('-new',optn)
ifig0=length(dir([fname_,'_SmoothedUnobserved*.fig']));
end
else
afig=dir([fname_,'_SmoothedUnobserved*.fig']);
for j=1:length(afig),
delete(afig(j).name);
end
end
if exist(options_.datafile)
instr = options_.datafile;
else
instr = ['load ' options_.datafile];
end
eval(instr);
if ~exist('T','var'),
temp = eval(deblank(options_.varobs(1,:)));
T=[1:length(temp)]';
clear temp;
end
fobs = options_.first_obs;
nv=length(varargin);
for ifig = 1:ceil(nv/9),
figure('Name','Smoothed unobserved variables'),
for j=1+9*(ifig-1):min(9*ifig,nv),
subplot(3,3,j-9*(ifig-1)),
i = strmatch(varargin{j},lgy_,'exact');
if ~isempty(texname)
if strmatch('\log',texname(i,:));
texname(i,1:end-1)=texname(i,2:end);
end
end
eval(['plot(T(fobs:fobs+options_.nobs-1),[oo_.SmoothedVariables.',varargin{j},'(1:end)+ys_(i)],''-k'')'])
if exist(varargin{j},'var')
hold on,
eval(['plot(T(fobs:fobs+options_.nobs-1),',varargin{j},'(fobs:fobs+options_.nobs-1),'':k'')'])
end
if options_.TeX,
title(deblank(texname(i,:)),'interpreter','tex')
else
title(varargin{j},'interpreter','none')
end
set(gca,'xlim',[T(fobs)-1 T(fobs+options_.nobs-1)+1])
end
saveas(gcf,[fname_,'_SmoothedUnobserved',int2str(ifig+ifig0)])
eval(['print -depsc2 ' fname_,'_SmoothedUnobserved',int2str(ifig+ifig0)]);
eval(['print -dpdf ' fname_,'_SmoothedUnobserved',int2str(ifig+ifig0)]);
if options_.nograph
close(gcf)
end
end

141
matlab/prior_draw_gsa.m
View File

@ -1,141 +0,0 @@
function pdraw = prior_draw_gsa(init,rdraw,cc)
% Build one draw from the prior distribution.
%
% INPUTS
% o init [integer] scalar equal to 1 (first call) or 0.
% o cc [double] two columns matrix (same as in
% metropolis.m), constraints over the
% parameter space (upper and lower bounds).
%
% OUTPUTS
% o pdraw [double] draw from the joint prior density.
%
% ALGORITHM
% ...
%
% SPECIAL REQUIREMENTS
% None.
%
%
% part of DYNARE, copyright S. Adjemian, M. Juillard (2006)
% Gnu Public License.
global M_ options_ estim_params_ bayestopt_
persistent fname npar bounds pshape pmean pstd a b p3 p4 condition
if init
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
MhDirectoryName = CheckPath('metropolis');
fname = [ MhDirectoryName '/' M_.fname];
pshape = bayestopt_.pshape;
pmean = bayestopt_.pmean;
pstd = bayestopt_.pstdev;
p1 = bayestopt_.p1;
p2 = bayestopt_.p2;
p3 = bayestopt_.p3;
p4 = bayestopt_.p4;
a = zeros(npar,1);
b = zeros(npar,1);
if nargin == 2
bounds = cc;
else
bounds = kron(ones(npar,1),[-Inf Inf]);
end
for i = 1:npar
switch pshape(i)
case 3% Gaussian prior
b(i) = pstd(i)^2/(pmean(i)-p3(i));
a(i) = (pmean(i)-p3(i))/b(i);
case 1% Beta prior
mu = (p1(i)-p3(i))/(p4(i)-p3(i));
stdd = p2(i)/(p4(i)-p3(i));
a(i) = (1-mu)*mu^2/stdd^2 - mu;
b(i) = a(i)*(1/mu - 1);
case 2;%Gamma prior
mu = p1(i)-p3(i);
b(i) = p2(i)^2/mu;
a(i) = mu/b(i);
case {5,4,6}
% Nothing to do here
%
% 4: Inverse gamma, type 1, prior
% p2(i) = nu
% p1(i) = s
% 6: Inverse gamma, type 2, prior
% p2(i) = nu
% p1(i) = s
% 5: Uniform prior
% p3(i) and p4(i) are used.
otherwise
disp('prior_draw :: Error!')
disp('Unknown prior shape.')
return
end
pdraw = zeros(npar,1);
end
condition = 1;
pdraw = zeros(npar,1);
return
end
for i = 1:npar
switch pshape(i)
case 5% Uniform prior.
pdraw(:,i) = rdraw(:,i)*(p4(i)-p3(i)) + p3(i);
case 3% Gaussian prior.
pdraw(:,i) = norminv(rdraw(:,i),pmean(i),pstd(i));
case 2% Gamma prior.
pdraw(:,i) = gaminv(rdraw(:,i),a(i),b(i))+p3(i);
case 1% Beta distribution (TODO: generalized beta distribution)
pdraw(:,i) = betainv(rdraw(:,i),a(i),b(i))*(p4(i)-p3(i))+p3(i);
case 4% INV-GAMMA1 distribution
while condition
tmp = sqrt(1/gamma_draw(p2(i)/2,1/p1(i),0));
if tmp >= bounds(i,1) && tmp <= bounds(i,2)
pdraw(i) = tmp;
break
end
end
case 6% INV-GAMMA2 distribution
while condition
tmp = 1/gamma_draw(p2(i)/2,1/p1(i),0);
if tmp >= bounds(i,1) && tmp <= bounds(i,2)
pdraw(i) = tmp;
break
end
end
otherwise
% Nothing to do here.
end
end
function g = gamma_draw(a,b,c)
% Bauwens, Lubrano & Richard (page 316)
if a >30
z = randn;
g = b*(z+sqrt(4*a-1))^2/4 + c;
else
condi = 1;
while condi
x = -1;
while x<0
u1 = rand;
y = tan(pi*u1);
x = y*sqrt(2*a-1)+a-1;
end
u2 = rand;
if log(u2) <= log(1+y^2)+(a-1)*log(x/(a-1))-y*sqrt(2*a-1);
break
end
end
g = x*b+c;
end

74
matlab/probam.m
View File

@ -1,74 +0,0 @@
function [paths_avg,res,bounds,res1]=probam(paths,prc,neigv)
%
% Input: paths(i,:) is the i^th path,
% 'prc' is a row vector containing percents (20=20%),
% 'neigv' is the third dimension of res(:,:,:), see below.
%
% Output: paths_avg is the average of paths (row vector).
% For each i, paths(i,:) can be expressed as
% paths_avg+sum_j{gamma(i,j)*eigvec(:,j)'}, where
% eigvec(:,j)' is the j^th eigenvector.
% For a given percent prc(l), let maki_lk denote the prc(l) percentile
% of gamma(:,k). Then res(l,:,k)=paths_avg+maki_lm*eigvec(:,m)', where m is
% the index of the eigenvector corresponding to the k^th largest eigenvalue.
% bounds is simple (see below).
% If m is the index of the eigenvector corresponding to the largest
% eigenvalue, and gamma(i,m) is the element of gamma(:,m) corresponding
% to the prc(l) percentile of gamma(:,m), then res1(l,:)=paths(i,:).
% Other details are in probam.pdf
[npaths,lpaths]=size(paths);
nprc=length(prc);
bounds=zeros(2,lpaths);
prcb=[prc(1) prc(nprc)];
for i=1:lpaths
[prctemp1,prctemp2]=prctil(paths(:,i),prcb);
bounds(:,i)=prctemp1';
end
w=cov(paths);
[eigvec,eigval]=eig(w);
eigvals=diag(eigval);
paths_avg=zeros(1,lpaths);
for i=1:npaths
paths_avg=paths_avg+paths(i,:);
end
paths_avg=paths_avg/npaths;
gamma=zeros(npaths,lpaths);
for i=1:npaths
gamma(i,:)=(eigvec'*(paths(i,:)-paths_avg)')';
end
[eigsrt,eigord]=sort(eigvals);
res=zeros(nprc,lpaths,neigv);
for i=1:neigv
[prctemp1,prctemp2]=prctil(gamma(:,eigord(lpaths-i+1)),prc);
res(:,:,i)=(eigvec(:,eigord(lpaths-i+1))*prctemp1)';
for j=1:nprc
res(j,:,i)=res(j,:,i)+paths_avg;
end
end
[prctemp1,prctemp2]=prctil(gamma(:,eigord(lpaths)),prc);
res1=paths(prctemp2(:),:);
function [prcx,prcxk]=prctil(prcy,prcz)
%prcx returns prcz percentiles of vector prcy, and prcxk
%is such that prcx=prcy(prcxk)
ly=length(prcy);
lz=length(prcz);
prcx=prcz;
prcxk=prcz;
[temp,ix]=sort(prcy); %temp=prcy(ix)
prcz=(ly/100)*prcz;
for i=1:lz
if prcz(i)<ly/2
prcz(i)=ceil(prcz(i));
else
prcz(i)=ly-ceil(ly-prcz(i))+1;
end
prcx(i)=temp(prcz(i)); %prcx=prcy(ix(prcz))
prcxk(i)=ix(prcz(i));
end

132
matlab/redform_map.m
View File

@ -1,132 +0,0 @@
function redform_map(anamendo, anamlagendo, anamexo, iload, pprior, ilog, threshold, ksstat, alpha2, dirname)
%function redform_map(namendo, namlagendo, namexo, iload, pprior, ilog, threshold, ksstat, alpha2, dirname)
% copyright Marco Ratto 2006
global M_ oo_ estim_params_ options_
pnames = M_.param_names(estim_params_.param_vals(:,1),:);
if nargin<4 | isempty(iload),
iload=0;
end
if nargin<5 | isempty(pprior),
pprior=1;
end
if nargin<6 | isempty(ilog),
ilog=0;
end
if nargin<7,
threshold=[];
end
if nargin<8,
ksstat=0.1;
end
if nargin<9,
alpha2=0.2;
end
if nargin<10,
dirname='';
end
if pprior
load([dirname,'\',M_.fname,'_prior']);
adir=[dirname '\redform_stab'];
else
load([dirname,'\',M_.fname,'_mc']);
adir=[dirname '\redform_mc'];
end
if isempty(dir(adir))
mkdir(adir)
end
adir0=pwd;
%cd(adir)
nspred=size(T,2)-M_.exo_nbr;
x0=lpmat(istable,:);
[kn, np]=size(x0);
nsok = length(find(M_.lead_lag_incidence(M_.maximum_lag,:)));
for j=1:size(anamendo,1)
namendo=deblank(anamendo(j,:));
iendo=strmatch(namendo,M_.endo_names(oo_.dr.order_var,:),'exact');
for jx=1:size(anamexo,1)
namexo=deblank(anamexo(jx,:));
iexo=strmatch(namexo,M_.exo_names,'exact');
if ~isempty(iexo),
%y0=squeeze(T(iendo,iexo+nspred,istable));
y0=squeeze(T(iendo,iexo+nspred,:));
clear lpmat lpmat0 egg iunstable yys
xdir0 = [adir,'\',namendo,'_vs_', namexo];
if isempty(threshold)
redform_private(x0, y0, iload, pnames, namendo, namexo, xdir0)
else
iy=find( (y0>threshold(1)) & (y0<threshold(2)));
iyc=find( (y0<=threshold(1)) | (y0>=threshold(2)));
xdir = [xdir0,'_cut'];
redform_private(x0(iy,:), y0(iy), iload, pnames, namendo, namexo, xdir)
delete([xdir, '\*cut*.*'])
[proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0);
indsmirnov = find(dproba>ksstat);
stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir);
stab_map_2(x0(iy,:),alpha2,'cut',xdir)
stab_map_2(x0(iyc,:),alpha2,'trim',xdir)
end
if ilog==-1,
yy=log(-y0);
xdir=[xdir0,'_minuslog'];
redform_private(x0, yy, iload, pnames, namendo, namexo, xdir)
elseif ilog==1,
yy=log(y0);
xdir=[xdir0,'_log'];
redform_private(x0, yy, iload, pnames, namendo, namexo, xdir)
end
end
end
for je=1:size(anamlagendo,1)
namlagendo=deblank(anamlagendo(je,:));
ilagendo=strmatch(namlagendo,M_.endo_names(oo_.dr.order_var(oo_.dr.nstatic+1:oo_.dr.nstatic+nsok),:),'exact');
if ~isempty(ilagendo),
%y0=squeeze(T(iendo,ilagendo,istable));
y0=squeeze(T(iendo,ilagendo,:));
xdir0 = [adir,'\',namendo,'_vs_', namlagendo];
if isempty(threshold)
redform_private(x0, y0, iload, pnames, namendo, namlagendo, xdir0)
else
iy=find( (y0>threshold(1)) & (y0<threshold(2)));
iyc=find( (y0<=threshold(1)) | (y0>=threshold(2)));
xdir = [xdir0,'_cut'];
redform_private(x0(iy,:), y0(iy), iload, pnames, namendo, namlagendo, xdir)
delete([xdir, '\*cut*.*'])
[proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0);
indsmirnov = find(dproba>ksstat);
stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir);
stab_map_2(x0(iy,:),alpha2,'cut',xdir)
stab_map_2(x0(iyc,:),alpha2,'trim',xdir)
end
if ilog==-1,
yy=log(-y0);
xdir=[xdir0,'_minuslog'];
redform_private(x0, yy, iload, pnames, namendo, namlagendo, xdir)
elseif ilog==1,
yy=log(y0);
xdir=[xdir0,'_log'];
redform_private(x0, yy, iload, pnames, namendo, namlagendo, xdir)
end
end
end
end
function redform_private(x0, y0, iload, pnames, namy, namx, xdir)
figure, hist(y0,30), title([namy,' vs. ', namx])
if isempty(dir(xdir))
mkdir(xdir)
end
saveas(gcf,[xdir,'\', namy,'_vs_', namx])
eval(['print -depsc2 ' xdir,'\', namy,'_vs_', namx]);
eval(['print -dpdf ' xdir,'\', namy,'_vs_', namx]);
fname=[xdir,'\map'];
if iload==0,
gsa_ = gsa_sdp_fn(y0, x0, -2, fname, pnames, 1);
else
gsa_ = gsa_sdp_fn(y0, x0, 0, fname, pnames, 1);
end

437
matlab/rividIRF.m
View File

@ -1,437 +0,0 @@
function rivid_irf(m, nroot, var_list_, iload)
% Metropolis-Hastings algorithm.
%
% INPUTS
% o type [char] string specifying the joint density of the
% deep parameters ('prior','posterior').
%
% OUTPUTS
% None (oo_ and plots).
%
%
% ALGORITHM
% None.
%
% SPECIAL REQUIREMENTS
% None.
%
%
% part of DYNARE, copyright S. Adjemian, M. Juillard (2006)
% Gnu Public License.
global options_ estim_params_ oo_ M_ dsge_prior_weight
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
offset = npar-np;
%load caz
%M_.Sigma_e=caz;
if nargin<4, iload=0; end,
%%
MaxNumberOfPlotPerFigure = 9;% The square root must be an integer!
nn = sqrt(MaxNumberOfPlotPerFigure);
DirectoryName = CheckPath('Output');
MhDirectoryName = CheckPath('GSA');
MAX_nirfs = ceil(options_.MaxNumberOfBytes/(options_.irf*length(oo_.steady_state)*M_.exo_nbr)/8)+50;
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
load([ MhDirectoryName '/' M_.fname '_prior'],'lpmat','istable')
x=lpmat(istable,:);
clear lpmat
NumberOfDraws=size(x,1);
B=NumberOfDraws; options_.B = B;
if ~iload
try delete([MhDirectoryName '\' M_.fname '_IRFs*']);
catch disp('No _IRFs files to be deleted!')
end
irun = 0;
irun2 = 0;
NumberOfIRFfiles = 1;
ifil2 = 1;
h = waitbar(0,'GSA (prior) IRFs...');
if B <= MAX_nruns
stock_param = zeros(B, npar);
else
stock_param = zeros(MAX_nruns, npar);
end
if B >= MAX_nirfs
stock_irf = zeros(options_.irf,M_.endo_nbr,M_.exo_nbr,MAX_nirfs);
else
stock_irf = zeros(options_.irf,M_.endo_nbr,M_.exo_nbr,B);
end
u=[zeros(10,1); 1; zeros(options_.irf-1,1)];
i0=[];
izp=[];
for b=1:B,
zp=0;
irun = irun+1;
irun2 = irun2+1;
deep = x(b,:);
stock_param(irun2,:) = deep;
set_parameters(deep);
%set_coefficients(deep); % only runs coefficients fixing the st. error of shocks
dr = resol(oo_.steady_state,0);
SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord) = M_.Sigma_e+1e-14*eye(M_.exo_nbr);
cs = 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:M_.endo_nbr,
jj=strmatch(deblank(M_.endo_names(j,:)),var_list_,'exact');
if jj,
yy=[zeros(10,1); y(j,:)'];
tf=squeeze(m(i,jj,:))';
[th,stats]=riv([yy u],tf,10);
if stats(3)<1-1.e-6,
disp('');
%[th,stats]=rivid([yy u],[[tf(1:2)-1 tf(3:4)]; [tf(1:2)+1 tf(3:4)]],1,10);
[a,bb,c,P,d]=getpar(th);
tf(1)=length(a)-1;
tf(3)=length(find(bb==0));
tf(2)=length(find(bb));
tf(4)=0;
end
sta(i,jj,b,:)=stats;
[a,bb,c,P,d]=getpar(th);
am{i,jj}(b,1:tf(1))=a(2:end);
bm{i,jj}(b,1:tf(2))=bb(tf(3)+1:end);
ar=roots(a);
for ir=1:tf(1),
if find(abs(roots(bb)-ar(ir))<1.e-5)
zp=1;
[th,stats]=riv([yy u],[tf(1)-1,tf(2)-1,tf(3:end)],10);
sta(i,jj,b,:)=stats;
[a,bb,c,P,d]=getpar(th);
am{i,jj}(b,:)=am{i,jj}(b,:)*0;
bm{i,jj}(b,:)=bm{i,jj}(b,:)*0;
am{i,jj}(b,1:tf(1)-1)=a(2:end);
bm{i,jj}(b,1:tf(2)-1)=bb(tf(3)+1:end);
end
end
end
if max(y(j,:)) - min(y(j,:)) > 1e-10
stock_irf(:,j,i,irun) = transpose(y(j,:));
end
end
end
end
if irun == MAX_nirfs | irun == B | b == B
if b == B
stock_irf = stock_irf(:,:,:,1:irun);
end
save([MhDirectoryName '/' M_.fname '_irf' int2str(NumberOfIRFfiles)],'stock_irf');
NumberOfIRFfiles = NumberOfIRFfiles+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)],'stock');
ifil2 = ifil2 + 1;
irun2 = 0;
end
if min(min(sta(:,:,b,3)))<(1-1.e-6), % | zp,
i0=[i0 b];
if zp,
izp=[izp b];
end
end,
waitbar(b/B,h);
end
NumberOfIRFfiles = NumberOfIRFfiles-1;
ifil2 = ifil2-1;
close(h);
ReshapeMatFiles('irf','gsa');
i2=find(~ismember([1:B],i0));
save([MhDirectoryName '/' M_.fname '_rivid'],'x','am','bm','i0','i2','izp','sta');
else
load([MhDirectoryName '/' M_.fname '_rivid'],'x','am','bm','i0','i2','izp','sta');
end
a0=[];
for i=1:M_.exo_nbr,
for j=1:size(var_list_,1)
a0=[a0 am{i,j}(i2,:)];
end
end
r0={};
r00=[];
a00=[];
%nroot=3;
izp=[];
vzp={};
h = waitbar(0,'Analysing roots of denominators...');
for b=1:length(i2),
roo=[];
zp=0;
wzp=[];
for i=1:M_.exo_nbr,
for j=1:size(var_list_,1),
tf=squeeze(m(i,j,:))';
bb=bm{i,j}(i2(b),:);
r0{i,j}(b,:)=sort(roots([1 am{i,j}(i2(b),:)]));
% bb=bm{i,j}(b,:);
% r0{i,j}(b,:)=sort(roots([1 am{i,j}(b,:)]));
for ir=1:tf(1),
if r0{i,j}(b,ir)~=0 & find(abs(roots(bb)-r0{i,j}(b,ir))<1.e-5)
zp=1;
wzp=[wzp; [i j]];
end
end
roo=[roo r0{i,j}(b,:)];
% if j>1 | i>1,
% f0=[];
% for l=1:m(i,j,1)
% f0(l)=find(abs(((r0{i,j}(b,:))-(r0{i,1}(b,l))))<=1.e-6);
% end
% r0{i,j}(b,:)=r0{i,j}(b,f0);
% end
end
end
roo=sort(roo);
crit=1;
rxx=[];
while length(rxx)~=nroot & zp==0 & crit>1.e-10,
crit=crit*0.5;
ix=find(abs(diff(sign(real(roo)).*abs(roo)))>crit);
ix=[1 ix+1];
rx=roo(ix);
rxx=[];
for jx=1:length(rx),
rxx=[rxx rx(jx)];
if find(imag(rx(jx))),
rxx=[rxx conj(rx(jx))];
end
end
end
if zp==0 & crit>1.e-10,
r00=[r00; rxx];
a00=[a00; poly(rxx)];
else
izp=[izp i2(b)];
% izp=[izp b];
vzp=[vzp; {wzp}];
r00=[r00; zeros(1,nroot)];
a00=[a00; poly(zeros(1,nroot))];
end
waitbar(b/length(i2),h)
% waitbar(b/1669,h)
end
close(h)
i1=find(~ismember([i2],izp));
save([MhDirectoryName '/' M_.fname '_rivid'],'r00','a00','-append');
if options_.opt_gsa.morris,
load([MhDirectoryName '/' M_.fname '_prior'],'lpmat');
lmat=[];
inx=[];
for j=1:options_.opt_gsa.morris_ntra,
z(j)=length(find(ismember([1+(estim_params_.np+1)*(j-1):(estim_params_.np+1)*j],izp)));
if z(j)==0,
inx=[inx [1+12*(j-1):12*j]];
lmat=[lmat; lpmat([1+12*(j-1):12*j],:)];
end
end,
i1=inx;
end
x=x(i1,:);
ya=a00(i1,2:end);
j0=size(ya,2);
for j=1:M_.exo_nbr,
j00=0;
yb=[];
yr=[];
y=am{j,1}(i1,:);
for l=1:size(y,2),
y0=y(:,l);
if max(y0)-min(y0)>1.e-10,
j00=j00+1;
yb=[yb y0];
end
end
for i=1:size(var_list_,1),
y=bm{j,i}(i1,:);
for l=1:size(y,2),
y0=y(:,l);
if max(y0)-min(y0)>1.e-10,
j00=j00+1;
yb=[yb y0];
j0=j0+1;
ya=[ya y0];
end
end
end
ybb(j)={yb};
if ~options_.opt_gsa.morris,
for jj=1:j00,
[dum, is]=sort(yb(:,jj));
yr(is,jj)=[1:length(i1)]'./length(i1);
end
ybr(j)={yr};
end
end
if options_.opt_gsa.morris==1,
nliv = options_.opt_gsa.morris_nliv;
for j=1:j0,
[SAmeas, SAMorris(:,:,j)] = Morris_Measure_Groups(estim_params_.np, x, ya(:,j),nliv);
end
SAM = squeeze(SAMorris(:,1,:));
for j=1:j0,
SAnorm(:,j)=SAM(:,j)./max(SAM(:,j));
end
figure, bar(SAnorm), title('All')
figure, bar(SAnorm'), title('All Relationships')
if exist('d0')
er=diag(d0);
%id=find((er./length(er))>0.01);
id=find(cumsum(er(end:-1:1))./j0<0.99);
id=length(id)+1;
jpc=0;
for j=j0:-1:(j0-id(1)+1),
jpc=jpc+1;
[SAmeas_r, SAMorris_r(:,:,jpc)] = Morris_Measure_Groups(estim_params_.np, x, ya*v0(:,j),nliv);
end
SAM = squeeze(SAMorris_r(:,1,:));
for j=1:jpc,
SAnorm_r(:,j)=SAM(:,j)./max(SAM(:,j));
end
figure, bar(SAnorm_r), title('All PCA')
for jj=1:27, yr(:,jj) = interp1(c.ya(:,jj),c.yr(:,jj),ya(:,jj)); end,
er=diag(d);
%id=find((er./length(er))>0.01);
id=find(cumsum(er(end:-1:1))./j0<0.99);
id=length(id)+1;
jpc=0;
for j=j0:-1:(j0-id(1)+1),
jpc=jpc+1;
[SAmeas_r, SAMorris_r(:,:,jpc)] = Morris_Measure_Groups(estim_params_.np, x, yr*v(:,j),nliv);
end
SAM = squeeze(SAMorris_r(:,1,:));
for j=1:jpc,
SAnorm_r(:,j)=SAM(:,j)./max(SAM(:,j));
end
figure, bar(SAnorm_r), title('All PCA rank')
end
for j=1:M_.exo_nbr,
SAMo=[];
yb = ybb{j};
for l=1:size(yb,2),
y0=yb(:,l);
[SAmeas, SAMo(:,:,l)] = Morris_Measure_Groups(estim_params_.np, x, y0,nliv);
end
SAM = squeeze(SAMo(:,1,:));
SAno=[];
for jj=1:size(yb,2),
SAno(:,jj)=SAM(:,jj)./max(SAM(:,jj));
end
figure, bar(SAno), title(M_.exo_names(j,:))
figure, bar(SAno'), title([M_.exo_names(j,:),' relationships' ])
end
elseif options_.opt_gsa.morris==0,
yr=[];
for j=1:j0,
[dum, is]=sort(ya(:,j));
yr(is,j)=[1:length(i1)]'./length(i1);
end
gsa_flag=0; %0; %-2 for GSA estimation
RividDir = CheckPath('GSA/rivid');
for j=1:j0,
gsa_(j) = gsa_sdp_fn(ya(:,j), x, gsa_flag, [RividDir,'/map',int2str(j)], M_.param_names);
end
% for j=1:j0,
% gsa_y(j) = gsa_sdp_fn(ya(:,j), ya(:,[1:j-1,j+1:j0]), 0, [RividDir,'/map_y',int2str(j)], M_.param_names);
% end
for j=1:j0, S(j,:)=gsa_(j).multivariate.si; end,
[v0,d0]=eig(corrcoef(ya));
save([MhDirectoryName '/' M_.fname '_rivid'],'v0','d0','-append');
ee=diag(d0);
%id=find((er./length(er))>0.01);
id=find(cumsum(ee(end:-1:1))./j0<0.99);
id=length(id)+1;
jpc=0;
for j=j0:-1:(j0-id(1)+1),
jpc=jpc+1;
gsa_pca(jpc) = gsa_sdp_fn(ya*v0(:,j), x, gsa_flag, [RividDir,'/map_pca',int2str(jpc)], M_.param_names);
end
for j=1:jpc, S_pca(j,:)=gsa_pca(j).multivariate.si; end,
for j=1:j0,
gsar_(j) = gsa_sdp_fn(yr(:,j), x, gsa_flag, [RividDir,'/map_rank',int2str(j)], M_.param_names);
end
for j=1:j0, Sr(j,:)=gsar_(j).multivariate.si; end,
[v,d]=eig(corrcoef(yr));
er=diag(d);
%id=find((er./length(er))>0.01);
id=find(cumsum(er(end:-1:1))./j0<0.99);
id=length(id)+1;
jpc=0;
for j=j0:-1:(j0-id(1)+1),
jpc=jpc+1;
gsar_pca(jpc) = gsa_sdp_fn(yr*v(:,j), x, gsa_flag, [RividDir,'/map_rank_pca',int2str(jpc)], M_.param_names);
end
for j=1:jpc, Sr_pca(j,:)=gsar_pca(j).multivariate.si; end,
fi=zeros(id,estim_params_.np);
for j=1:id,
yup=find(Sr_pca(j,:)>0.05);
if ~isempty(yup),
fi(j,yup)=ones(1,length(yup));
end
end
% for j=1:j0-id,
% gsar_pca_0(j) = gsa_sdp_fn(yr*v(:,j), x, gsa_flag, [RividDir,'/rivid/map_rank_pca_0',int2str(j)], M_.param_names);
% end
% for j=1:j0-id, Sr_pca_0(j,:)=gsar_pca_0(j).multivariate.si; end,
% yr0=zeros(length(i1),1);
% for j=1:10, %j0-id,
% yr0 = yr0+yr*v(:,j);
% end
% gsar_pca_00 = gsa_sdp_fn(yr0, x, gsa_flag, [RividDir,'/map_rank_pca_00'], M_.param_names);
% Sr_pca_00=gsar_pca_00.multivariate.si;
% y0 = getIRFRuns(1,1,'gsa');
% prc=[2.5 16 84 97.5];
% [paths_avg,res,bounds,res1]=probam(y0',prc,2);
% figure,
% plot([1:options_.irf],paths_avg,'k-',[1:options_.irf],bounds(1,:),'k:', [1:options_.irf], bounds(2,:),'k:');
save([MhDirectoryName '/' M_.fname '_rivid'],'v','d','ya','yr','-append');
jpc=0;
for j=j0:-1:(j0-id(1)+1),
jpc=jpc+1;
y0=yr*v(:,j);
rbf_1(jpc)=rbf(x, y0, 1);
disp(['Output ',num2str(jpc),' out of ',num2str(id(1)),' done!'])
end
jpc=0;
for j=j0:-1:(j0-id(1)+1),
jpc=jpc+1;
y0=yr*v(:,j);
ivec=find(gsar_pca(jpc).multivariate.si>0.01);
y0=y0-sum(gsar_pca(jpc).multivariate.f(:,ivec),2);
rbf_2(jpc)=rbf(x(:,ivec), y0, 2);
disp(['Output ',num2str(jpc),' out of ',num2str(id(1)),' done!'])
end
save([MhDirectoryName '/' M_.fname '_rivid'],'rbf_1','rbf_2','-append');
end

257
matlab/sdr.m
View File

@ -1,257 +0,0 @@
function [fit,fitse,par,parse,zs,pars,parses,rsq,nvr,y0]=...
sdp0(y,z,x,TVP,nvr,opts,P0,x0,nvr0,tab)
% SDP Non-parametric state dependent regression modelling (backfitting)
%
% [fit,fitse,par,parse,xs,pars,parses,rsq,nvre,y0]=...
% sdp(y,z,x,TVP,nvr,opts,P0,x0,nvr0,tab)
% 1 2 3 4 5 6 7 8 9
%
% y: Time series (*)
% z: Regressors (*)
% x: States on which the parameters depend, column for each regressor (z)
% TVP: Model type for each TVP (0-RW, 1-IRW) (0)
% nvr: NVR hyper-parameters (0)
% if an nvr is a negative integer (-n) it is
% optimised over the first n backfit steps
% opts: estimation options [iter con meth sm ALG plotopt], use -1 for defaults
% iter: number of backfitting iterations (10)
% con: backfitting convergence threshold (0.001)
% meth: Optimisation estimation method (0)
% 0: Maximum Likelihood
% Integer:Sum of squares of the #-step-ahead forecasting errors
% sm: Smoothing on (1-default) or off (0-saves memory)
% ALG: Smoothing algorithm P (0) or Q (1-default)
% plotopt: Plot results during estimation on (1) or off (0-default)
% P0: Initial P matrix diagonal (1e5)
% x0: Initial state vector (0)
% nvr0: Initial NVRs specified by user (0.0001)
% tab: controls progress display, 0 - no display, 2 - continuous display
%
% fit: Model fit
% fitse: Standard error of the fit
% par: Parameter estimates
% parse: Standard errors of parameters
% xs: Sorted states
% pars: Sorted parameter estimates
% parses: Sorted standard errors of parameters
% rsq: R squared value
% nvre: Estimated NVRs
% y0: Interpolated data
%
% Example: sdr(y,[u1 u2],[x1 x2],[0 1],[-2 -1])
% regression type model y = c1(x1)*u1 + c2(x1)*u2, with an RW model for
% c1 where the dependent state is x1 and an IRW model for c2 where the
% dependent state is x2; the NVR for the first SDP (c1) is optimised at
% the first iteration and at the first two iterations for the second SDP
%
% See also FCAST, STAND, DLR, DHR
%
% Copyright (c) 2004 by CRES, Lancaster University, United Kingdom
% Authors : Peter Young, Wlodek Tych, Diego Pedregal, James Taylor
% Additional author: Paul McKenna
% Captain Tbx. functions used: dlr01, ksirw2c, dlropt01, dlrfun01
% Revision history :
% 28/01/2005, WT, conversion into SDR solver, code cosmetics
% 24/09/2003, PM, corrected for x0 input with IRW
% 23/09/2003, PM, corrected for x0 input
% 05/03/2003, JT, interpolated data output argument
% 22/10/2002, PM, adapted for convergence with missing values in output
% 10/07/2002, JT, default NVR changed to 1 in line nvr(j)=min(nvr2,1)
% 30/05/2002, JT, included in toolbox, P0 now passed to DLROPT as a vector
% 18/04/2002, PM, nvr input made more flexible for optimisation
% 17/04/2002, PM, optimisation set to ALG=0 to avoid problems with for Matlab 4
% 17/04/2002, PM, nvr input padded with end values if too short
% 15/04/2002, PM, opts input argument created to make the function call more compact
% Argument sizes
Ndat=size(z,1);nr=size(z,2);
% ******* Default argument settings
% NVR0
if nargin<10,tab=[]; end
if nargin<9;nvr0=[];end
if isempty(nvr0);nvr0=0.0001;end
if length(nvr0)==1 & nr~=1;nvr0=nvr0*ones(1,nr);end
lnvr0=length(nvr0);nvr0=nvr0(:)';nvr0=[nvr0 ones(1,nr-lnvr0)*nvr0(lnvr0)];
% x0: initial conditions (state)
if nargin<8,x0=[];end
% P0: initial conditions on P-matrix
if nargin<7,P0=1e5;end
if isempty(P0);P0=1e5;end
% options
if nargin<6,opts=[];end
if isempty(opts),opts=[10 0.001 0 1 1 0];end
opts0=[10 0.001 0 1 1 0];
if length(opts)~=6;
opts=[opts opts0(length(opts)+1:6)];
end
if isempty(tab), tab=2; end
% ******* Options setting
opts(find(opts<0))=opts0(find(opts<0));
bf_it=opts(1);
con=opts(2);
sm=opts(4);
ALG=opts(5);
plotopt=opts(6);
if opts(3)==0;
meth='ml';
else
meth=['f' int2str(opts(3))];
end
% ********
% Default NVR settings
if nargin<5,nvr=0;end
if isempty(nvr);nvr=0;end
if length(nvr)==1 & nr~=1;nvr=nvr*ones(1,nr);end
lnvr=length(nvr);nvr=nvr(:)';nvr=[nvr ones(1,nr-lnvr)*nvr(lnvr)];
% Default TVP settings (type of RW model)
if nargin<4,TVP=0;end
if isempty(TVP);TVP=0;end
if length(TVP)==1 & nr~=1;TVP=TVP*ones(1,nr);end
lTVP=length(TVP);TVP=TVP(:)';TVP=[TVP ones(1,nr-lTVP)*TVP(lTVP)];
if ~isempty(P0),lP0=length(P0);P0=P0(:)';P0=[P0 ones(1,sum(TVP)+nr-lP0)*P0(lP0)];end
% Default settings for state and regressors
if nargin <3, x=z;end
if isempty(x);x=z;end
%marking the state locations that relate to nvrs (not alpha)
ipp=[];
for i=1:length(TVP);ipp=[ipp 1 zeros(1,TVP(i))];end
ipp=find(ipp);
optim=zeros(size(nvr));
par=zeros(Ndat,nr);
parse=zeros(Ndat,nr);fit=zeros(Ndat,1);
ccv=zeros(bf_it,1); % convergence criteria vector
for i=1:length(nvr);
if nvr(i)<0;
optim(i)=-nvr(i);
end
end
% Estimate time invariant parameter estimates
par0=(z\y)';
par=par0(ones(length(y),1),:);
% Optional graphics output
if plotopt;
clf;
for i=1:nr;
subplot(nr,1,i);
%plot(z(:,i),par(:,i),'b-');drawnow;hold on
ll=plot(x(:,i),par(:,i),'b.');set(ll,'markersize',1);drawnow;hold on
end
end
% Optimisation options - variable recycling
opts=foptions;
opts(2)=0.01;%termination tolerance
i=0;
for i=1:bf_it;
for j=1:nr;
par0=par;
% Generate partial residual
npar=[1:nr];npar=[npar(npar<j) npar(npar>j)];
if nr==1;
yc=y;
else
if nr>2;
yc=y-sum(par0(:,npar)'.*z(:,npar)')';
else
yc=y-(par0(:,npar).*z(:,npar));
end
end
% Sort data
[ys,I]=sort(x(:,j)); %sort w.r.t. j-th regressor
ysc=yc(I); %partial residual, sorted and corrected
xr=z(:,j); %regressor
xs=xr(I); %regressor, sorted
%optimise
if i<=optim(j);
if i==1;
[nvr2,alpha2,opts2,parse2]=dlropt01(ysc,xs,TVP(j),meth,-2,1,nvr0(j),1,opts,0,tab,P0(ipp(j):ipp(j)+TVP(j)));% ALG=0 for Matlab4
else
[nvr2,alpha2,opts2,parse2]=dlropt01(ysc,xs,TVP(j),meth,-2,1,nvr(j),1,opts,0,tab,P0(ipp(j):ipp(j)+TVP(j)));% ALG=0 for Matlab4
end
nvr(j)=min(nvr2,1); % default=1
end
% Estimate
if i==1
if isempty(x0);
[fit,fitse,par1,parse1,comp,e]=dlr01(ysc,xs,TVP(j),nvr(j),1,P0(ipp(j):ipp(j)+TVP(j)),x0,1);
else
[fit,fitse,par1,parse1,comp,e]=dlr01(ysc,xs,TVP(j),nvr(j),1,P0(ipp(j):ipp(j)+TVP(j)),[x0(j);zeros(TVP(j),1)],1);
end
else
[fit,fitse,par1,parse1,comp,e]=dlr01(ysc,xs,TVP(j),nvr(j),1,P0(ipp(j):ipp(j)+TVP(j)),[par(I(1),j);zeros(TVP(j),1)],1);
end
par(I,j)=par1;
parse(I,j)=parse1;
if plotopt;
subplot(nr,1,j);
if i>1;eval(['set(l' int2str(j) ',''color'',[0.5 0.5 0.5]);']);end
eval(['l' int2str(j) '=plot(x(:,j),par(:,j),''r.'');']);eval(['set(l' int2str(j) ',''markersize'',1);']);drawnow
end
end
if i<=max(optim);
disp(['Optimised NVRs (after ' int2str(i) ' of ' int2str(bf_it) ' backfitting iterations)'])
disp(nvr)
if all(optim<=i) & j==nr;
disp('Retaining these values for subsequent iterations')
end
end
NotMiss=find(~isnan(y) & ~sum(isnan(par)')' & ~sum(isnan(z)')' );% Flag "not missing" values
if nr>1;
disp([i 1-(cov(y(NotMiss)-sum(par(NotMiss,:)'.*z(NotMiss,:)')')/cov(y(NotMiss)))])
else
disp([i 1-(cov(y(NotMiss)-(par(NotMiss,:).*z(NotMiss,:)))/cov(y(NotMiss)))])
end
if nr>1;
ccv(i)=1-(cov(y(NotMiss)-sum(par(NotMiss,:)'.*z(NotMiss,:)')')/cov(y(NotMiss)));
else
ccv(i)=1-(cov(y(NotMiss)-(par(NotMiss,:).*z(NotMiss,:)))/cov(y(NotMiss)));
end
if i>3;
if abs(ccv(i)-ccv(i-1))<con | abs(ccv(i)-ccv(i-2))<con;
disp(['Converged on backfit iteration number ' int2str(i)]);break
else
if i==bf_it;
disp('Maximum number of backfit iterations reached, convergence criterion not achieved')
end
end
end
end
[zs,zI]=sort(x);
pars=zeros(size(par));
parses=zeros(size(parse));
for i=1:nr;
pars(:,i)=par(zI(:,i),i);
parses(:,i)=parse(zI(:,i),i);
end
if nr>1;
fit=sum(par'.*z')';
else
fit=par.*z;
end
nvre=nvr;
rsq=1-(cov(y(NotMiss)-fit(NotMiss))/cov(y(NotMiss)));
% interpolated data
ii=find(isnan(y));
y0=y;
y0(ii)=fit(ii);
% end of m-file

BIN
matlab/sdr.p
View File

Binary file not shown.

14
matlab/set_coefficients.m
View File

@ -1,14 +0,0 @@
function set_coefficients(xparam1)
global estim_params_ M_
nvx = estim_params_.nvx;
ncx = estim_params_.ncx;
nvn = estim_params_.nvn;
ncn = estim_params_.ncn;
np = estim_params_.np;
if np
offset = nvx+ncx+nvn+ncn;
M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end);
end

58
matlab/smirnov.m
View File

@ -1,58 +0,0 @@
function [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
% Smirnov test for 2 distributions
% [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
%
% Copyright (C) 2005 Marco Ratto
%
if nargin<3
alpha = 0.05;
end
if nargin<4,
iflag=0;
end
% empirical cdfs.
xmix= [x1;x2];
bin = [-inf ; sort(xmix) ; inf];
ncount1 = histc (x1 , bin);
ncount2 = histc (x2 , bin);
cum1 = cumsum(ncount1)./sum(ncount1);
cum1 = cum1(1:end-1);
cum2 = cumsum(ncount2)./sum(ncount2);
cum2 = cum2(1:end-1);
n1= length(x1);
n2= length(x2);
n = n1*n2 /(n1+n2);
% Compute the d(n1,n2) statistics.
if iflag==0,
d = max(abs(cum1 - cum2));
elseif iflag==-1
d = max(cum2 - cum1);
elseif iflag==1
d = max(cum1 - cum2);
end
%
% Compute P-value check H0 hypothesis
%
lam = max((sqrt(n) + 0.12 + 0.11/sqrt(n)) * d , 0);
if iflag == 0
j = [1:101]';
prob = 2 * sum((-1).^(j-1).*exp(-2*lam*lam*j.^2));
prob=max(prob,0);
prob=min(1,prob);
else
prob = exp(-2*lam*lam);
end
H = (alpha >= prob);

606
matlab/stab_map_.m
View File

@ -1,606 +0,0 @@
function x0 = stab_map_(Nsam, fload, ksstat, alpha2, prepSA, pprior, ilptau, OutputDirectoryName)
%
% function x0 = stab_map_(Nsam, fload, alpha2, prepSA, pprior)
%
% Mapping of stability regions in the prior ranges applying
% Monte Carlo filtering techniques.
%
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models
% I. Mapping stability, MIMEO, 2005.
%
% INPUTS
% Nsam = MC sample size
% fload = 0 to run new MC; 1 to load prevoiusly generated analysis
% alpha2 = significance level for bivariate sensitivity analysis
% [abs(corrcoef) > alpha2]
% prepSA = 1: save transition matrices for mapping reduced form
% = 0: no transition matrix saved (default)
% pprior = 1: sample from prior ranges (default): sample saved in
% _prior.mat file
% = 0: sample from posterior ranges: sample saved in
% _mc.mat file
% OUTPUT:
% x0: one parameter vector for which the model is stable.
%
% GRAPHS
% 1) Pdf's of marginal distributions under the stability (dotted
% lines) and unstability (solid lines) regions
% 2) Cumulative distributions of:
% - stable subset (dotted lines)
% - unacceptable subset (solid lines)
% 3) Bivariate plots of significant correlation patterns
% ( abs(corrcoef) > alpha2) under the stable and unacceptable subsets
%
% USES lptauSEQ,
% stab_map_1, stab_map_2
%
% Copyright (C) 2005 Marco Ratto
% THIS PROGRAM WAS WRITTEN FOR MATLAB BY
% Marco Ratto,
% Unit of Econometrics and Statistics AF
% (http://www.jrc.cec.eu.int/uasa/),
% IPSC, Joint Research Centre
% The European Commission,
% TP 361, 21020 ISPRA(VA), ITALY
% marco.ratto@jrc.it
%
% ALL COPIES MUST BE PROVIDED FREE OF CHARGE AND MUST INCLUDE THIS COPYRIGHT
% NOTICE.
%
%global bayestopt_ estim_params_ dr_ options_ ys_ fname_
global bayestopt_ estim_params_ options_ oo_ M_
opt_gsa=options_.opt_gsa;
nliv = opt_gsa.morris_nliv;
ntra = opt_gsa.morris_ntra;
dr_ = oo_.dr;
if isfield(dr_,'ghx'),
ys_ = oo_.dr.ys;
nspred = size(dr_.ghx,2);
nboth = dr_.nboth;
nfwrd = dr_.nfwrd;
end
fname_ = M_.fname;
nshock = estim_params_.nvx;
nshock = nshock + estim_params_.nvn;
nshock = nshock + estim_params_.ncx;
nshock = nshock + estim_params_.ncn;
lpmat0=[];
if nargin==0,
Nsam=2000; %2^13; %256;
end
if nargin<2,
fload=0;
end
if nargin<3,
ksstat=0.1;
end
if nargin<4,
alpha2=0.3;
end
if nargin<5,
prepSA=0;
end
if nargin<6,
pprior=1;
end
if nargin<7,
ilptau=1;
end
if nargin<8,
OutputDirectoryName='';
end
options_.periods=0;
options_.nomoments=1;
options_.irf=0;
options_.noprint=1;
options_.simul=0;
if fload==0 | nargin<2 | isempty(fload),
if prepSA
T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam/2);
end
if isfield(dr_,'ghx'),
egg=zeros(length(dr_.eigval),Nsam);
end
yys=zeros(length(dr_.ys),Nsam);
if opt_gsa.morris
if opt_gsa.morris == 1
[lpmat, OutFact] = Sampling_Function_2(nliv, estim_params_.np, ntra, ones(estim_params_.np, 1), zeros(estim_params_.np,1), []);
lpmat = lpmat.*(nliv-1)/nliv+1/nliv/2;
Nsam=size(lpmat,1);
elseif opt_gsa.morris==2
lpmat = prep_ide(Nsam,estim_params_.np,5);
Nsam=size(lpmat,1);
end
else
if estim_params_.np<52 & ilptau>0,
[lpmat] = lptauSEQ(Nsam,estim_params_.np); % lptau
if estim_params_.np>30 | ilptau==2, % scrambled lptau
for j=1:estim_params_.np,
lpmat(:,j)=lpmat(randperm(Nsam),j);
end
end
else ilptau==0
%[lpmat] = rand(Nsam,estim_params_.np);
for j=1:estim_params_.np,
lpmat(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
end
end
end
prior_draw_gsa(1);
if pprior,
for j=1:nshock,
lpmat0(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
%lpmat0(:,j)=lpmat0(:,j).*(bayestopt_.ub(j)-bayestopt_.lb(j))+bayestopt_.lb(j);
end
%for j=1:estim_params_.np,
% lpmat(:,j)=lpmat(:,j).*(bayestopt_.ub(j+nshock)-bayestopt_.lb(j+nshock))+bayestopt_.lb(j+nshock);
%end
xx=prior_draw_gsa(0,[lpmat0 lpmat]);
lpmat0=xx(:,1:nshock);
lpmat=xx(:,nshock+1:end);
clear xx;
else
% for j=1:nshock,
% xparam1(j) = oo_.posterior_mode.shocks_std.(bayestopt_.name{j});
% sd(j) = oo_.posterior_std.shocks_std.(bayestopt_.name{j});
% lpmat0(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
% lb = max(bayestopt_.lb(j), xparam1(j)-2*sd(j));
% ub1=xparam1(j)+(xparam1(j) - lb); % define symmetric range around the mode!
% ub = min(bayestopt_.ub(j),ub1);
% if ub<ub1,
% lb=xparam1(j)-(ub-xparam1(j)); % define symmetric range around the mode!
% end
% lpmat0(:,j) = lpmat0(:,j).*(ub-lb)+lb;
% end
% %
% for j=1:estim_params_.np,
% xparam1(j+nshock) = oo_.posterior_mode.parameters.(bayestopt_.name{j+nshock});
% sd(j+nshock) = oo_.posterior_std.parameters.(bayestopt_.name{j+nshock});
% lb = max(bayestopt_.lb(j+nshock),xparam1(j+nshock)-2*sd(j+nshock));
% ub1=xparam1(j+nshock)+(xparam1(j+nshock) - lb); % define symmetric range around the mode!
% ub = min(bayestopt_.ub(j+nshock),ub1);
% if ub<ub1,
% lb=xparam1(j+nshock)-(ub-xparam1(j+nshock)); % define symmetric range around the mode!
% end
% %ub = min(bayestopt_.ub(j+nshock),xparam1(j+nshock)+2*sd(j+nshock));
% if estim_params_.np>30 & estim_params_.np<52
% lpmat(:,j) = lpmat(randperm(Nsam),j).*(ub-lb)+lb;
% else
% lpmat(:,j) = lpmat(:,j).*(ub-lb)+lb;
% end
% end
%load([fname_,'_mode'])
eval(['load ' options_.mode_file ';']');
d = chol(inv(hh));
lp=randn(Nsam*2,nshock+estim_params_.np)*d+kron(ones(Nsam*2,1),xparam1');
for j=1:Nsam*2,
lnprior(j) = any(lp(j,:)'<=bayestopt_.lb | lp(j,:)'>=bayestopt_.ub);
end
ireal=[1:2*Nsam];
ireal=ireal(find(lnprior==0));
lp=lp(ireal,:);
Nsam=min(Nsam, length(ireal));
lpmat0=lp(1:Nsam,1:nshock);
lpmat=lp(1:Nsam,nshock+1:end);
clear lp lnprior ireal;
end
%
h = waitbar(0,'Please wait...');
istable=[1:Nsam];
jstab=0;
iunstable=[1:Nsam];
iindeterm=zeros(1,Nsam);
iwrong=zeros(1,Nsam);
for j=1:Nsam,
M_.params(estim_params_.param_vals(:,1)) = lpmat(j,:)';
stoch_simul([]);
dr_ = oo_.dr;
if isfield(dr_,'ghx'),
egg(:,j) = sort(dr_.eigval);
iunstable(j)=0;
if prepSA
jstab=jstab+1;
T(:,:,jstab) = [dr_.ghx dr_.ghu];
end
if ~exist('nspred'),
nspred = size(dr_.ghx,2);
nboth = dr_.nboth;
nfwrd = dr_.nfwrd;
end
else
istable(j)=0;
if isfield(dr_,'eigval')
egg(:,j) = sort(dr_.eigval);
if exist('nspred')
if any(isnan(egg(1:nspred,j)))
iwrong(j)=j;
else
if (nboth | nfwrd) & abs(egg(nspred+1,j))<=options_.qz_criterium,
iindeterm(j)=j;
end
end
end
else
if exist('egg'),
egg(:,j)=ones(size(egg,1),1).*1.1;
end
iwrong(j)=j;
end
end
ys_=real(dr_.ys);
yys(:,j) = ys_;
ys_=yys(:,1);
waitbar(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
end
close(h)
if prepSA,
T=T(:,:,1:jstab);
end
istable=istable(find(istable)); % stable params
iunstable=iunstable(find(iunstable)); % unstable params
iindeterm=iindeterm(find(iindeterm)); % indeterminacy
iwrong=iwrong(find(iwrong)); % dynare could not find solution
% % map stable samples
% istable=[1:Nsam];
% for j=1:Nsam,
% if any(isnan(egg(1:nspred,j)))
% istable(j)=0;
% else
% if abs(egg(nspred,j))>=options_.qz_criterium; %(1-(options_.qz_criterium-1)); %1-1.e-5;
% istable(j)=0;
% %elseif (dr_.nboth | dr_.nfwrd) & abs(egg(nspred+1,j))<=options_.qz_criterium; %1+1.e-5;
% elseif (nboth | nfwrd) & abs(egg(nspred+1,j))<=options_.qz_criterium; %1+1.e-5;
% istable(j)=0;
% end
% end
% end
% istable=istable(find(istable)); % stable params
%
% % map unstable samples
% iunstable=[1:Nsam];
% for j=1:Nsam,
% %if abs(egg(dr_.npred+1,j))>1+1.e-5 & abs(egg(dr_.npred,j))<1-1.e-5;
% %if (dr_.nboth | dr_.nfwrd),
% if ~any(isnan(egg(1:5,j)))
% if (nboth | nfwrd),
% if abs(egg(nspred+1,j))>options_.qz_criterium & abs(egg(nspred,j))<options_.qz_criterium; %(1-(options_.qz_criterium-1));
% iunstable(j)=0;
% end
% else
% if abs(egg(nspred,j))<options_.qz_criterium; %(1-(options_.qz_criterium-1));
% iunstable(j)=0;
% end
% end
% end
% end
% iunstable=iunstable(find(iunstable)); % unstable params
if pprior,
if ~prepSA
save([OutputDirectoryName '\' fname_ '_prior'],'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','nspred','nboth','nfwrd')
else
save([OutputDirectoryName '\' fname_ '_prior'],'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','T','nspred','nboth','nfwrd')
end
else
if ~prepSA
save([OutputDirectoryName '\' fname_ '_mc'], ...
'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','nspred','nboth','nfwrd')
else
save([OutputDirectoryName '\' fname_ '_mc'], ...
'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','T','nspred','nboth','nfwrd')
end
end
else
if pprior,
filetoload=[OutputDirectoryName '\' fname_ '_prior'];
else
filetoload=[OutputDirectoryName '\' fname_ '_mc'];
end
load(filetoload,'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','nspred','nboth','nfwrd')
Nsam = size(lpmat,1);
if prepSA & isempty(strmatch('T',who('-file', filetoload),'exact')),
h = waitbar(0,'Please wait...');
options_.periods=0;
options_.nomoments=1;
options_.irf=0;
options_.noprint=1;
stoch_simul([]);
T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),length(istable));
ntrans=length(istable);
for j=1:ntrans,
M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(j),:)';
stoch_simul([]);
dr_ = oo_.dr;
T(:,:,j) = [dr_.ghx dr_.ghu];
if ~exist('nspred')
nspred = size(dr_.ghx,2);
nboth = dr_.nboth;
nfwrd = dr_.nfwrd;
end
ys_=real(dr_.ys);
yys(:,j) = ys_;
ys_=yys(:,1);
waitbar(j/ntrans,h,['MC iteration ',int2str(j),'/',int2str(ntrans)])
end
close(h)
save(filetoload,'T','-append')
else
load(filetoload,'T')
end
end
if pprior
aname='prior_stab';
auname='prior_unacceptable';
aunstname='prior_unstable';
aindname='prior_indeterm';
asname='prior_stable';
else
aname='mc_stab';
auname='mc_unacceptable';
aunstname='mc_unstable';
aindname='mc_indeterm';
asname='mc_stable';
end
delete([OutputDirectoryName,'\',fname_,'_',aname,'_*.*']);
%delete([OutputDirectoryName,'\',fname_,'_',aname,'_SA_*.*']);
delete([OutputDirectoryName,'\',fname_,'_',asname,'_corr_*.*']);
delete([OutputDirectoryName,'\',fname_,'_',auname,'_corr_*.*']);
delete([OutputDirectoryName,'\',fname_,'_',aunstname,'_corr_*.*']);
delete([OutputDirectoryName,'\',fname_,'_',aindname,'_corr_*.*']);
if length(iunstable)>0 & length(iunstable)<Nsam,
disp([num2str(length(istable)/Nsam*100),'\% of the prior support is stable.'])
disp([num2str( (length(iunstable)-length(iwrong)-length(iindeterm) )/Nsam*100),'\% of the prior support is unstable.'])
if ~isempty(iindeterm),
disp([num2str(length(iindeterm)/Nsam*100),'\% of the prior support gives indeterminacy.'])
end
if ~isempty(iwrong),
disp(' ');
disp(['For ',num2str(length(iwrong)/Nsam*100),'\% of the prior support dynare could not find a solution.'])
end
% Blanchard Kahn
[proba, dproba] = stab_map_1(lpmat, istable, iunstable, aname,0);
indstab=find(dproba>ksstat);
disp('The following parameters mostly drive acceptable behaviour')
disp(M_.param_names(estim_params_.param_vals(indstab,1),:))
stab_map_1(lpmat, istable, iunstable, aname, 1, indstab, OutputDirectoryName);
if ~isempty(iindeterm),
ixun=iunstable(find(~ismember(iunstable,[iindeterm,iwrong])));
[proba, dproba] = stab_map_1(lpmat, [1:Nsam], iindeterm, [aname, '_indet'],0);
indindet=find(dproba>ksstat);
disp('The following parameters mostly drive indeterminacy')
disp(M_.param_names(estim_params_.param_vals(indindet,1),:))
stab_map_1(lpmat, [1:Nsam], iindeterm, [aname, '_indet'], 1, indindet, OutputDirectoryName);
if ~isempty(ixun),
[proba, dproba] = stab_map_1(lpmat, [1:Nsam], ixun, [aname, '_unst'],0);
indunst=find(dproba>ksstat);
disp('The following parameters mostly drive instability')
disp(M_.param_names(estim_params_.param_vals(indunst,1),:))
stab_map_1(lpmat, [1:Nsam], ixun, [aname, '_unst'], 1, indunst, OutputDirectoryName);
end
end
disp(' ')
disp('Starting bivariate analysis:')
c0=corrcoef(lpmat(istable,:));
c00=tril(c0,-1);
stab_map_2(lpmat(istable,:),alpha2, asname, OutputDirectoryName);
stab_map_2(lpmat(iunstable,:),alpha2, auname, OutputDirectoryName);
if ~isempty(iindeterm),
stab_map_2(lpmat(iindeterm,:),alpha2, aindname, OutputDirectoryName);
if ~isempty(ixun),
stab_map_2(lpmat(ixun,:),alpha2, aunstname, OutputDirectoryName);
end
end
x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
x0 = [x0; lpmat(istable(1),:)'];
if istable(end)~=Nsam
M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(1),:)';
stoch_simul([]);
end
else
if length(iunstable)==0,
disp('All parameter values in the specified ranges are stable!')
x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
x0 = [x0; lpmat(istable(1),:)'];
else
disp('All parameter values in the specified ranges are not acceptable!')
x0=[];
end
end
if opt_gsa.redform,
[nr,nc,nn]=size(T);
j0=0;
for j=1:nr,
for i=1:nc,
y0=squeeze(T(j,i,:));
if max(y0)-min(y0)>1.e-10,
j0=j0+1;
y1=ones(size(lpmat,1),1)*NaN;
y1(istable,1)=y0;
yt(:,j0)=y1;
end
end
end
end
if opt_gsa.morris==1 & opt_gsa.redform,
OutputDir = CheckPath('GSA');
for j=1:j0,
[SAmeas, SAMorris(:,:,j)] = Morris_Measure_Groups(estim_params_.np, lpmat, yt(:,j),nliv);
end
SAM = squeeze(SAMorris(:,1,:));
for j=1:j0
SAnorm(:,j)=SAM(:,j)./max(SAM(:,j));
irex(j)=length(find(SAnorm(:,j)>0.01));
end
[dum, irel]=sort(irex);
figure, bar(SAnorm(:,irel))
set(gca,'xtick',[1:estim_params_.np])
set(gca,'xlim',[0.5 estim_params_.np+0.5])
title('Elementary effects parameters')
saveas(gcf,[OutputDir,'\',fname_,'_morris_par'])
eval(['print -depsc2 ' OutputDir '\' fname_ '_morris_par']);
eval(['print -dpdf ' OutputDir '\' fname_ '_morris_par']);
figure, bar(SAnorm(:,irel)')
set(gca,'xtick',[1:j0])
set(gca,'xlim',[0.5 j0+0.5])
title('Elementary effects relationships')
saveas(gcf,[OutputDir,'\',fname_,'_morris_redform'])
eval(['print -depsc2 ' OutputDir '\' fname_ '_morris_redform']);
eval(['print -dpdf ' OutputDir '\' fname_ '_morris_redform']);
elseif opt_gsa.morris==2 & opt_gsa.redform,
np=estim_params_.np;
na=(4*np+1)*opt_gsa.Nsam;
for j=1:j0,
[idex(j,:), yd(j,:)] = spop_ide(lpmat, yt(:,j), opt_gsa.Nsam, 5-1);
end
iok=find(~isnan(yt(1:opt_gsa.Nsam,1)));
yr=NaN*ones(size(lpmat,1),j0);
for j=1:j0,
ys(j,:)=yd(j,:)./max(yd(j,:));
[dum, is]=sort(yt(iok,j));
yr(iok(is),j)=[1:length(iok)]'./length(iok);
yr(istable(length(iok)+1:end),j) = interp1(yt(iok,j),yr(iok,j),yt(istable(length(iok)+1:end),j),'','extrap');
ineg=find(yr(:,j)<0);
if any(ineg),
[dum, is]=sort(yr(ineg,j));
yr(ineg(is),j)=-[length(ineg):-1:1]./length(iok);
end
[idex_r(j,:), yd_r(j,:)] = spop_ide(lpmat, yr(:,j), opt_gsa.Nsam, 5-1);
ys_r(j,:)=yd_r(j,:)./max(yd_r(j,:));
end,
figure, bar((idex.*ys)./opt_gsa.Nsam), title('Relationships')
figure, bar((idex.*ys)'./opt_gsa.Nsam), title('Parameters')
figure, bar((idex_r.*ys_r)./opt_gsa.Nsam), title('Relationships rank')
figure, bar((idex_r.*ys_r)'./opt_gsa.Nsam), title('Parameters rank')
[v0,d0]=eig(corrcoef(yt(iok,:)));
ee=diag(d0);
ee=ee([end:-1:1])./j0;
i0=length(find(ee>0.01));
v0=v0(:,[end:-1:1]);
for j=1:i0,
[idex_pc(j,:), yd_pc(j,:)] = spop_ide(lpmat, yt*v0(:,j), opt_gsa.Nsam, 5-1);
end
for j=1:i0,
ys_pc(j,:)=yd_pc(j,:)./max(yd_pc(j,:));
end,
figure, bar((idex_pc.*ys_pc)./opt_gsa.Nsam), title('Relationships PCA')
figure, bar((idex_pc.*ys_pc)'./opt_gsa.Nsam), title('Parameters PCA')
[vr,dr]=eig(corrcoef(yr(iok,:)));
er=diag(dr);
er=er([end:-1:1])./j0;
ir0=length(find(er>0.01));
vr=vr(:,[end:-1:1]);
for j=1:ir0,
[idex_pcr(j,:), yd_pcr(j,:)] = spop_ide(lpmat, yr*vr(:,j), opt_gsa.Nsam, 5-1);
end
for j=1:ir0,
ys_pcr(j,:)=yd_pcr(j,:)./max(yd_pcr(j,:));
end,
figure, bar((idex_pcr.*ys_pcr)./opt_gsa.Nsam), title('Relationships rank PCA')
figure, bar((idex_pcr.*ys_pcr)'./opt_gsa.Nsam), title('Parameters rank PCA')
elseif opt_gsa.redform
yr=[];
for j=1:j0,
[dum, is]=sort(yt(istable,j));
yr(is,j)=[1:length(istable)]'./length(istable);
end
gsa_flag=0; %-2 for GSA estimation
RividDir = CheckPath('GSA/rivid');
for j=1:j0,
gsa_(j) = gsa_sdp_fn(yt(istable,j), lpmat(istable,:), ...
bayestopt_.pshape, [bayestopt_.p1 bayestopt_.p2 bayestopt_.p3 bayestopt_.p4], ...
gsa_flag, [RividDir,'/map_T',int2str(j)], M_.param_names);
gsa_r(j) = gsa_sdp_fn(yr(:,j), lpmat(istable,:), ...
bayestopt_.pshape, [bayestopt_.p1 bayestopt_.p2 bayestopt_.p3 bayestopt_.p4], ...
gsa_flag, [RividDir,'/map_T_rank',int2str(j)], M_.param_names);
S(j,:)=gsa_(j).multivariate.si;
Sr(j,:)=gsa_r(j).multivariate.si;
end
figure, bar(S')
set(gca,'xtick',[1:estim_params_.np])
set(gca,'xlim',[0.5 estim_params_.np+0.5])
title('Main effects parameters')
saveas(gcf,[RividDir,'\',fname_,'_SA_par'])
eval(['print -depsc2 ' RividDir '\' fname_ '_SA_par']);
eval(['print -dpdf ' RividDir '\' fname_ '_SA_par']);
figure, bar(S)
set(gca,'xtick',[1:j0])
set(gca,'xlim',[0.5 j0+0.5])
title('Main effects relationships')
saveas(gcf,[RividDir,'\',fname_,'_SA_redform'])
eval(['print -depsc2 ' RividDir '\' fname_ '_SA_redform']);
eval(['print -dpdf ' RividDir '\' fname_ '_SA_redform']);
figure, bar(Sr')
set(gca,'xtick',[1:estim_params_.np])
set(gca,'xlim',[0.5 estim_params_.np+0.5])
title('Main effects parameters rank')
saveas(gcf,[RividDir,'\',fname_,'_SA_par'])
eval(['print -depsc2 ' RividDir '\' fname_ '_SA_par_rank']);
eval(['print -dpdf ' RividDir '\' fname_ '_SA_par']);
figure, bar(Sr)
set(gca,'xtick',[1:j0])
set(gca,'xlim',[0.5 j0+0.5])
title('Main effects relationships rank')
saveas(gcf,[RividDir,'\',fname_,'_SA_redform_rank'])
eval(['print -depsc2 ' RividDir '\' fname_ '_SA_redform_rank']);
eval(['print -dpdf ' RividDir '\' fname_ '_SA_redform_rank']);
[v0,d0]=eig(corrcoef(yt(istable,:)));
ee=diag(d0);
ee=ee(end:-1:1);
v0=v0(:,end:-1:1);
%id=find((er./length(er))>0.01);
id=find(cumsum(ee)./j0<0.99);
jpc=length(id)+1;
gsa_flag=-2;
for j=1:jpc,
gsa_pca(j) = gsa_sdp_fn(yt(istable,:)*v0(:,j), lpmat(istable,:), ...
bayestopt_.pshape, [bayestopt_.p1 bayestopt_.p2 bayestopt_.p3 bayestopt_.p4], ...
gsa_flag, [RividDir,'/map_T_pca',int2str(j)], M_.param_names);
S_pca(j,:)=gsa_pca(j).multivariate.si;
end
figure, bar(S_pca')
set(gca,'xtick',[1:estim_params_.np])
set(gca,'xlim',[0.5 estim_params_.np+0.5])
title('Main effects parameters PCA')
saveas(gcf,[RividDir,'\',fname_,'_SA_par_PCA'])
eval(['print -depsc2 ' RividDir '\' fname_ '_SA_par_PCA']);
eval(['print -dpdf ' RividDir '\' fname_ '_SA_par_PCA']);
figure, bar(S_pca)
set(gca,'xtick',[1:j0])
set(gca,'xlim',[0.5 j0+0.5])
title('Main effects relationships PCA')
saveas(gcf,[RividDir,'\',fname_,'_SA_redform_PCA'])
eval(['print -depsc2 ' RividDir '\' fname_ '_SA_redform_PCA']);
eval(['print -dpdf ' RividDir '\' fname_ '_SA_redform_PCA']);
end

71
matlab/stab_map_1.m
View File

@ -1,71 +0,0 @@
function [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, iplot, ipar, dirname)
%function [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, iplot, ipar, dirname)
%
% lpmat = Monte Carlo matrix
% ibehaviour = index of behavioural runs
% inonbehaviour = index of non-behavioural runs
% aname = label of the analysis
% iplot = 1 plot cumulative distributions (default)
% iplot = 0 no plots
% ipar = index array of parameters to plot
% dirname = (OPTIONAL) path of the directory where to save
% (default: current directory)
%
% Plots: dotted lines for BEHAVIOURAL
% solid lines for NON BEHAVIOURAL
% USES smirnov
global estim_params_ bayestopt_ M_ options_
if nargin<5,
iplot=1;
end
fname_ = M_.fname;
if nargin<7,
dirname='';;
end
nshock = estim_params_.nvx;
nshock = nshock + estim_params_.nvn;
nshock = nshock + estim_params_.ncx;
nshock = nshock + estim_params_.ncn;
npar=size(lpmat,2);
ishock= npar>estim_params_.np;
if nargin<6 | isempty(ipar),
ipar=[1:npar];
end
nparplot=length(ipar);
% Smirnov test for Blanchard;
for j=1:npar,
[H,P,KSSTAT] = smirnov(lpmat(ibehaviour,j),lpmat(inonbehaviour,j));
proba(j)=P;
dproba(j)=KSSTAT;
end
if iplot
lpmat=lpmat(:,ipar);
ftit=bayestopt_.name(ipar+nshock*(1-ishock));
for i=1:ceil(nparplot/12),
figure('name',aname),
for j=1+12*(i-1):min(nparplot,12*i),
subplot(3,4,j-12*(i-1))
if ~isempty(ibehaviour),
h=cumplot(lpmat(ibehaviour,j));
set(h,'color',[0 0 0], 'linestyle',':')
end
hold on,
if ~isempty(inonbehaviour),
h=cumplot(lpmat(inonbehaviour,j));
set(h,'color',[0 0 0])
end
title([ftit{j},'. D-stat ', num2str(dproba(ipar(j)),2)],'interpreter','none')
end
saveas(gcf,[dirname,'\',fname_,'_',aname,'_SA_',int2str(i)])
eval(['print -depsc2 ' dirname '\' fname_ '_' aname '_SA_' int2str(i)]);
eval(['print -dpdf ' dirname '\' fname_ '_' aname '_SA_' int2str(i)]);
if options_.nograph, close(gcf), end
end
end

92
matlab/stab_map_2.m
View File

@ -1,92 +0,0 @@
%function stab_map_2(x,alpha2,istab,fnam)
function stab_map_2(x,alpha2,fnam, dirname)
% function stab_map_2(x,alpha2,fnam)
%
% Copyright (C) 2005 Marco Ratto
% THIS PROGRAM WAS WRITTEN FOR MATLAB BY
% Marco Ratto,
% Unit of Econometrics and Statistics AF
% (http://www.jrc.cec.eu.int/uasa/),
% IPSC, Joint Research Centre
% The European Commission,
% TP 361, 21020 ISPRA(VA), ITALY
% marco.ratto@jrc.it
%
%global bayestopt_ estim_params_ dr_ options_ ys_ fname_
global bayestopt_ estim_params_ options_ oo_ M_
npar=size(x,2);
ishock= npar>estim_params_.np;
if nargin<3,
fnam='';
end
if nargin<4,
dirname='';
end
ys_ = oo_.dr.ys;
dr_ = oo_.dr;
fname_ = M_.fname;
nshock = estim_params_.nvx;
nshock = nshock + estim_params_.nvn;
nshock = nshock + estim_params_.ncx;
nshock = nshock + estim_params_.ncn;
c0=corrcoef(x);
c00=tril(c0,-1);
fig_nam_=[fname_,'_',fnam,'_corr_'];
ifig=0;
j2=0;
if ishock==0
npar=estim_params_.np;
else
npar=estim_params_.np+nshock;
end
for j=1:npar,
i2=find(abs(c00(:,j))>alpha2);
if length(i2)>0,
for jx=1:length(i2),
j2=j2+1;
if mod(j2,12)==1,
ifig=ifig+1;
figure('name',['Correlations in the ',fnam,' sample ', num2str(ifig)]),
end
subplot(3,4,j2-(ifig-1)*12)
% bar(c0(i2,j)),
% set(gca,'xticklabel',bayestopt_.name(i2)),
% set(gca,'xtick',[1:length(i2)])
%plot(stock_par(ixx(nfilt+1:end,i),j),stock_par(ixx(nfilt+1:end,i),i2(jx)),'.k')
%hold on,
plot(x(:,j),x(:,i2(jx)),'.')
% xlabel(deblank(estim_params_.param_names(j,:)),'interpreter','none'),
% ylabel(deblank(estim_params_.param_names(i2(jx),:)),'interpreter','none'),
if ishock,
xlabel(bayestopt_.name{j},'interpreter','none'),
ylabel(bayestopt_.name{i2(jx)},'interpreter','none'),
else
xlabel(bayestopt_.name{j+nshock},'interpreter','none'),
ylabel(bayestopt_.name{i2(jx)+nshock},'interpreter','none'),
end
title(['cc = ',num2str(c0(i2(jx),j))])
if (mod(j2,12)==0) & j2>0,
saveas(gcf,[dirname,'\',fig_nam_,int2str(ifig)])
eval(['print -depsc2 ' dirname '\' fig_nam_ int2str(ifig)]);
eval(['print -dpdf ' dirname '\' fig_nam_ int2str(ifig)]);
if options_.nograph, close(gcf), end
end
end
end
if (j==(npar)) & j2>0,
saveas(gcf,[dirname,'\',fig_nam_,int2str(ifig)])
eval(['print -depsc2 ' dirname '\' fig_nam_ int2str(ifig)]);
eval(['print -dpdf ' dirname '\' fig_nam_ int2str(ifig)]);
if options_.nograph, close(gcf), end
end
end
if ifig==0,
disp(['No correlation term >', num2str(alpha2),' found for ',fnam])
end
%close all

152
matlab/stab_map_cyc.m
View File

@ -1,152 +0,0 @@
function x0 = stab_map_cyc(Period, alpha2, alpha, pprior)
%
% function x0 = stab_map_(Period, alpha2, alpha)
%
% Mapping of target cyclicity region in the prior ranges applying
% Monte Carlo filtering techniques
%
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models
% I. Mapping stability, MIMEO, 2005.
%
% INPUTS
% Period = range of periodicity for the targete cycles [Tmin Tmax]
% alpha2 = significance level for bivariate sensitivity analysis
% [abs(corrcoef) > alpha2]
% alpha = magnitude if cyclical eigenvalue w.r.t. max eigenvalue
% (default=0.9)
%
% OUTPUT:
% x0: one parameter vector for which the model provides targeted cycle.
%
% GRAPHS
% 1) Pdf's of marginal distributions under the ciclicity (dotted
% lines) and non-cyclicity (solid lines) regions
% 2) Cumulative distributions of:
% - cyclic subset (dotted lines)
% - uncyclic subset (solid lines)
% 3) Bivariate plots of significant correlation patterns
% ( abs(corrcoef) > alpha2) under the cyclic/uncyclic subsets
%
%
% Copyright (C) 2005 Marco Ratto
% THIS PROGRAM WAS WRITTEN FOR MATLAB BY
% Marco Ratto,
% Unit of Econometrics and Statistics AF
% (http://www.jrc.cec.eu.int/uasa/),
% IPSC, Joint Research Centre
% The European Commission,
% TP 361, 21020 ISPRA(VA), ITALY
% marco.ratto@jrc.it
%
% ALL COPIES MUST BE PROVIDED FREE OF CHARGE AND MUST INCLUDE THIS COPYRIGHT
% NOTICE.
%
% USES: stab_map_1
% stab_map_2
%global bayestopt_ estim_params_ dr_ options_ ys_ fname_
global bayestopt_ estim_params_ options_ oo_ M_
dr_ = oo_.dr;
if isfield(dr_,'ghx'),
ys_ = oo_.dr.ys;
nspred = size(dr_.ghx,2);
nboth = dr_.nboth;
nfwrd = dr_.nfwrd;
end
fname_ = M_.fname;
nshock = estim_params_.nvx;
nshock = nshock + estim_params_.nvn;
nshock = nshock + estim_params_.ncx;
nshock = nshock + estim_params_.ncn;
if nargin==0,
Nsam=2000; %2^13; %256;
end
if nargin<4,
pprior=1;
end
if nargin<3,
alpha=0.9;
end
if nargin<2,
alpha2=0.3;
end
if pprior,
load([fname_,'_prior'])
delete([fname_,'_cyc*.*']);
delete([fname_,'_cyclic*.*']);
else
load([fname_,'_mc'])
delete([fname_,'_cyc_mc*.*']);
delete([fname_,'_cyclic_mc*.*']);
end
Nsam = size(lpmat,1);
if length(istable)>0,
thex=istable;
for j=1:length(istable),
PerX{j}=[];
MaxEig = max(abs(egg(1:nspred,istable(j))));
ic = find(imag(egg(1:nspred,istable(j))));
%i=find( abs(egg( ic ,istable(j)) )>0.9*MaxEig); %only consider complex dominant eigenvalues
i=find( abs(egg( ic ,istable(j)) )>alpha*MaxEig); %only consider complex dominant eigenvalues
if ~isempty(i),
i=i(1:2:end);
thedum=[];
for ii=1:length(i),
idum = ic( i(ii) );
thedum(ii)=2*pi/abs(angle(egg(idum,istable(j))));
end
if thedum<Period(1) | thedum>Period(2),
thex(j)=0;
else
PerX{j}=[PerX{j} thedum(find(thedum>Period(1) & thedum<Period(2)))];
end
% [dum, icx]=max(thedum);
% icy(j) = ic( i(icx) );
% PerX(j)=max(thedum);
% if PerX(j)<Period(1) | PerX(j)>Period(2),
% thex(j)=0;
% end
else
thex(j)=0;
end
end
inocyc=istable(find(thex==0)); % non-cyclic params
icyc=thex(find(thex)); % cyclic params
if ~isempty(icyc)
if pprior,
cycnam='cyc';
cycnam2='cyclic';
else
cycnam='cyc_mc';
cycnam2='cyclic_mc';
end
proba = stab_map_1(lpmat, icyc, inocyc, cycnam);
stab_map_2(lpmat(icyc,:),alpha2, cycnam2);
x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
x0 = [x0; lpmat(icyc(1),:)'];
else
disp('None of the parameter values gives the targeted cyclicity!')
x0=[];
end
%stab_map_2(lpmat(iunstable,:),alpha2, 0);
else
if length(iunstable)==0,
disp('All parameter values in the specified ranges are stable!')
else
disp('All parameter values in the specified ranges are unstable!')
end
end

25
matlab/stand_.m
View File

@ -1,25 +0,0 @@
function [y, meany, stdy] = stand_(x)
% STAND_ Standardise a matrix by columns
%
% [x,my,sy]=stand(y)
%
% y: Time series (column matrix)
%
% x: standardised equivalent of y
% my: Vector of mean values for each column of y
% sy: Vector of standard deviations for each column of y
%
% Copyright (c) 2006 by JRC, European Commission, United Kingdom
% Author : Marco Ratto
if nargin==0,
return;
end
meany=mean(x);
stdy=std(x);
for j=1:size(x,2);
y(:,j)=(x(:,j)-meany(j))./stdy(j);
end
% end of m-file

BIN
matlab/stand_.p
View File

Binary file not shown.

98
matlab/vardec.m
View File

@ -1,98 +0,0 @@
function [Vi, Vk, Vkk] =vardec(k, varargin);
global oo_ M_
lgx_orig_ord_=M_.exo_names_orig_ord;
Sigma_e_=M_.Sigma_e;
lgy_=M_.endo_names;
lgx_=M_.exo_names;
fname_=M_.fname;
dr_=oo_.dr;
if isempty(k), k=40; end,
nvar = length(varargin);
if nvar == 0
nvar = length(dr_.order_var);
ivar = transpose(1:nvar);
var_list = lgy_;
else
ivar=zeros(nvar,1);
for i=1:nvar
if i>1,
var_list=str2mat(var_list,varargin{i});
else
var_list(i,:)=varargin{i};
end
i_tmp = strmatch(var_list(i,:),lgy_,'exact');
if isempty(i_tmp)
error (['One of the variable specified does not exist']) ;
else
ivar(i) = i_tmp;
end
end
end
[T,B,SteadyState,info] = dynare_resolve;
[dum, idum]=sort(dr_.order_var);
%BB=B(idum,:);
S0=diag(diag(Sigma_e_));
V=0.*T;
Vi=zeros(size(lgy_,1), size(lgx_,1), k);
Vii=zeros(size(T,1), size(T,2),size(lgx_,1));
%T=dr_.ghx;
%B=dr_.ghu;
for i=1:k,
V=T*V*T'+B*Sigma_e_*B';
V0=diag(V);
V0=V0(idum);
Vdum = V(idum, idum);
Vk(:,i)=V0(ivar);
Vkk(:,:,i)=Vdum(ivar, ivar);
iv{i}=find(V0(ivar)>1.e-12);
iv0=find(V0<=1.e-12);
V0(iv0)=ones(length(iv0),1);
for j=1:length(Sigma_e_),
Vii(:,:,j)= T*Vii(:,:,j)*T' + B(:,j)*S0(j,j)*B(:,j)';
V0i=diag(Vii(:,:,j));
Vi(:,j,i)=V0i(idum)./V0;
end
end
BQ=B*Sigma_e_*B';
Vinf = lyapunov_symm(T,BQ);
V0=diag(Vinf);
V0=V0(idum);
ivinf=find(V0(ivar)>1.e-12);
iv0=find(V0<=1.e-12);
V0(iv0)=ones(length(iv0),1);
for j=1:length(Sigma_e_),
BQ=B(:,j)*S0(j,j)*B(:,j)';
V00 = lyapunov_symm(T,BQ);
V00 = diag(V00);
Vinfi(:,j) = V00(idum)./V0;
%Vinfi(:,j) = inv(eye(length(T))-kron(T,T))*BQ(:)./V0;
end
iprint=[1, 4, 8, 16, 40, 100, 200, 400];
lh = size(deblank(lgy_(ivar,:)),2)+2;
headers = lgx_;
%headers(lgx_orig_ord_,:) = headers;
headers = strvcat(' ',headers);
for j=1:length(iprint),
ip=iprint(j);
if k>=ip,
title=['VARIANCE DECOMPOSITION (',int2str(ip),'-step ahead, in percent)'];
table(title,headers,deblank(lgy_(ivar(iv{ip}),:)),100*Vi(ivar(iv{ip}),:,ip), ...
lh,8,2);
end
end
title=['VARIANCE DECOMPOSITION (inf-step ahead, in percent)'];
table(title,headers,deblank(lgy_(ivar(ivinf),:)),100*Vinfi(ivar(ivinf),:), ...
lh,8,2);
save([fname_,'_vardec'],'Vi', 'Vinfi')

92
matlab/vecm0.m
View File

@ -1,92 +0,0 @@
function [bigVAR, yy1, yy2]=vecm0(xx, xvec, beta, nex, nk)
global opts_vecm_
a=ones(1,size(xx,2));
b=a*0;
b(nex)=ones(length(nex),1);
c=a-b;
nendo=find(c);
nvar=size(xx,2);
nvec=size(xvec,2);
for j=1:nvar;
if opts_vecm_.thmean
x0(j) = opts_vecm_.y0(j);
else
x0(j) = mean(xx(:,j));
end
xx(:,j)=xx(:,j)-x0(j);
end
if nvec,
x0vec=mean(xvec);
xvec=xvec-ones(length(xvec),1)*mean(xvec);
else
x0vec=[];
end
bigVAR=zeros(nvar,nvar+nvec);
% exogenous sub-block
xxx = [xx(1:end-1,nex)];
yyy = xx(2:end,nex) ;
bigVAR(nex,nex) = (inv(xxx'*xxx)*xxx'*yyy)';
if nvec,
[j0, i0]=find(beta(:,nex));
if ~isempty(j0),
if length(i0)==1,
xxx = [xx(1:end-1,nex) xvec(1:end-1,j0)];
bigVAR(nex(i0),[nex nvar+j0])= (inv(xxx'*xxx)*sum(xxx.*kron(yyy(:,(i0)),ones(1,size(xxx,2))))')';
elseif ~find(diff(i0))
xxx = [xx(1:end-1,nex) xvec(1:end-1,j0)];
bigVAR(nex(i0),[nex nvar+j0])= (inv(xxx'*xxx)*sum(xxx.*kron(yyy(:,(i0)),ones(1,size(xxx,2))))')';
else
jj=[1; find(diff(i0))+1];
ii=i0(jj);
for j=1:length(ii),
jj=find(i0==ii);
xxx = [xx(1:end-1,nex) xvec(1:end-1,j0(jj))];
bigVAR(nex(ii(jj)),[nex nvar+j0(jj)])=(inv(xxx'*xxx)*sum(xxx.*kron(yyy((ii(jj)),j),ones(1,size(xxx,2))))')';
end
end
end
end
% full VECM
if nvec,
xxx=[xx(1:end-1,:) xvec(1:end-1,:)];
else
xxx=[xx(1:end-1,:)];
end
yyy=xx(2:end,nendo);
for j=1:size(yyy,2),
bigVAR(nendo(j),:) = (inv(xxx'*xxx)*sum(xxx.*kron(yyy(:,j),ones(1,size(xxx,2))))')';
end
if nvec,
bigVAR = [bigVAR; ...
beta*bigVAR(:,1:nvar), (eye(size(beta,1))+beta*bigVAR(:,nvar+1:end) )]; % equilibrium correcton raws
end
%xxx=[xx(1:end-1,:) xvec(1:end-1,:)];
yy=(bigVAR*xxx')'; % 1-step ahead prediction
if nvec,
yyy=[xx(2:end,:) xvec(2:end,:)]; % data
else
yyy=xx(2:end,:); % data
end
xxx=[xx xvec];
x00=[x0 x0vec];
yy1=ones(size(xxx,1)+nk,size(xxx,2),nk);
for j=1:length(x00),
yy1(:,j,:)=yy1(:,j,:).*x00(j);
end
for j=1:nk,
yy1(j+1:size(xxx,1)+j,:, j)=(bigVAR^j*xxx')' + ones(size(xxx,1),1)*[x0 x0vec];
end
for j=1:nk,
yy2(j,:)=(bigVAR^j*xxx(end,:)')'+[x0 x0vec];
end

131
matlab/vecm2.m
View File

@ -1,131 +0,0 @@
function [bigVAR, yy1, yy2]=vecm2(xx, xvec, beta, nex, nk, nlag)
a=ones(1,size(xx,2));
b=a*0;
b(nex)=ones(length(nex),1);
c=a-b;
nendo=find(c);
nvar=size(xx,2);
nvec=size(xvec,2);
% for j=1:nvar;
% x0(j) = mean(xx(2:end,j));
% xx(:,j)=xx(:,j)-x0(j);
% end
%
% x0vec=mean(xvec(2:end,:));
% xvec=xvec-ones(length(xvec),1)*mean(xvec(2:end,:));
%
% exogenous sub-block
ind_ex_=1;
xxx=ones(size(xx,1)-nlag,1);
for j=1:nlag,
xxx = [xxx xx(j:end-nlag+j-1,nex)];
ind_ex_ = [ind_ex_ 1+nex+nvar*(j-1)];
end
yyy = xx(nlag+1:end,nex) ;
bigVAR=zeros(nvar,nvar*nlag+nvec+1);
for j=1:length(nex),
%bigrho(j,:) = (inv(xxx'*xxx)*sum(xxx.*kron(yyy(:,j),ones(1,size(xxx,2))))')';
bigVAR(nex(j),ind_ex_) = (inv(xxx'*xxx)*sum(xxx.*kron(yyy(:,j),ones(1,size(xxx,2))))')';
end
if nvec,
[j0, i0]=find(beta(:,nex));
if ~isempty(j0),
if length(i0)==1,
xxxx = [xxx xvec(1:end-nlag,j0)];
bigVAR(nex(i0),[ind_ex_ nvar*nlag+1+j0])= (inv(xxxx'*xxxx)*sum(xxxx.*kron(yyy(:,i0),ones(1,size(xxxx,2))))')';
elseif ~find(diff(i0))
xxxx = [xxx xvec(1:end-nlag,j0)];
bigVAR(nex(i0),[ind_ex_ nvar*nlag+1+j0])= (inv(xxxx'*xxxx)*sum(xxxx.*kron(yyy(:,i0),ones(1,size(xxxx,2))))')';
else
jj=[1; find(diff(i0))+1];
ii=i0(jj);
for j=1:length(ii),
jj=find(i0==ii);
%xxx = [xx(1:end-1,nex) xvec(1:end-1,j0(jj))];
xxxx = [xxx xvec(1:end-nlag,j0(jj))];
%bigVAR(i0,[nex nvar+j0(jj)])=(inv(xxx'*xxx)*sum(xxx.*kron(yyy(:,j),ones(1,size(xxx,2))))')';
bigVAR(nex(ii(jj)),[ind_ex_ nvar*nlag+1+j0(jj)])= (inv(xxxx'*xxxx)*sum(xxxx.*kron(yyy(:,ii(jj)),ones(1,size(xxxx,2))))')';
end
end
end
end
%bigrho = (inv(xxx'*xxx)*xxx'*yyy)';
% full VECM
xxx=ones(size(xx,1)-nlag,1);
for j=1:nlag,
xxx=[xxx xx(j:end-nlag+j-1,:)];
end
if nvec,
xxx=[xxx xvec(1:end-nlag,:)];
end
ind_vec_ = [nvar*nlag+2:nvar*nlag+nvec+1];
yyy=xx(nlag+1:end,nendo);
for j=1:length(nendo),
bigVAR(nendo(j),:) = (inv(xxx'*xxx)*sum(xxx.*kron(yyy(:,j),ones(1,size(xxx,2))))')';
end
if nvec,
dumVECM=beta*bigVAR(:,1);
for j=1:nlag,
dumVECM = [dumVECM beta*bigVAR(:,2+nvar*(j-1):nvar*j+1)];
end
dumVECM=[dumVECM eye(size(beta,1))+beta*bigVAR(:,ind_vec_) ];
if nlag>1,
bigVAR = [bigVAR; ...
[zeros(nvar*(nlag-1), 1) eye(nvar*(nlag-1)) zeros(nvar*(nlag-1), nvec+nvar)]; ... % more than 1 lag rows
dumVECM]; % equilibrium correcton raws
else
bigVAR = [bigVAR; ...
dumVECM]; % equilibrium correcton raws
end
end
yy=(bigVAR*xxx')'; % 1-step ahead prediction
yy=yy(:,1:(nvar+nvec));
if nvec,
yyy=[xx(nlag+1:end,:) xvec(nlag+1:end,:)]; % data
else
yyy=[xx(nlag+1:end,:)]; % data
end
dum0 = inv(eye(nlag*nvar+nvec)-bigVAR(:,2:end))*bigVAR(:,1);
x00=dum0([[1:nvar] ind_vec_-1]);
xxx=[];
for j=1:nlag,
%xxx=[xxx [zeros(j-1,nvar); ...
% xx(j:end,:)] ];
xxx=[xxx [ones(j-1,1)*x00(1:nvar)'; ...
xx(j:end,:)] ];
end
xxx=[xxx xvec];
yy1=ones(size(xxx,1)+nk,nvar+nvec,nk);
for j=1:length(x00),
yy1(:,j,:)=yy1(:,j,:).*x00(j);
end
for j=1:nk,
dum=xxx;
for i=1:j,
dum=[ones(size(xx,1),1) dum];
dum=(bigVAR*dum')';
end
yy1(j+1:size(xxx,1)+j,:, j)=dum(:,[[1:nvar] ind_vec_-1]);
end
for j=1:nk,
dum=xxx(end,:);
for i=1:j,
dum=[1 dum];
dum=(bigVAR*dum(end,:)')';
end
yy2(j,:)=dum([[1:nvar] ind_vec_-1]);
end
if max(sort(abs(eig(bigVAR(:,2:end)))))>1,
disp('WARNING: the VECM is unstable!!!')
end

122
matlab/vecm_prep.m
View File

@ -1,122 +0,0 @@
function [rmseVECM, rmseKVECM] = vecm_prep(varargin);
global options_ bayestopt_ M_ oo_
global opts_vecm_
ys_ = oo_.steady_state;
lgy_ = M_.endo_names;
fname_ = M_.fname;
beta=opts_vecm_.beta;
gflag=opts_vecm_.gflag;
if exist(options_.datafile)
instr = options_.datafile;
else
instr = ['load ' options_.datafile];
end
eval(instr);
if ~exist('T'), T=[1:options_.nobs]'; end
istart = options_.presample+1;
fobs = options_.first_obs;
nobs= options_.nobs;
if any(gflag) & opts_vecm_.lflag,
istart=istart-1;
fobs=fobs+1;
nobs=nobs-1;
end
xx=[];
x0=[];
for j=1:length(varargin),
eval(['dum=',varargin{j},';'])
x0=[x0 dum];
if gflag(j),
dum=[0; diff(dum)];
end
xx=[xx dum];
end
%xx = [E_GC E_GE E_GEX E_GI E_GIM E_GL E_GY E_WPHI E_INOM ...
% E_PHI E_PHIC E_PHIM E_PHIX E_GTR E_INOMW E_PHIW E_GYW];
xx=xx(fobs:fobs+nobs-1,:);
xvec=[];
for j=1:size(beta,1),
xvec = [xvec x0*beta(j,:)'];
end
%xvec = [ E_LCSN E_LISN E_LEXY E_LIMY E_LPCP E_LPMP E_LPXP E_L E_LWRY E_LTRYN E_LYWY];
if ~isempty(xvec)
xvec=xvec(fobs:fobs+nobs-1,:);
end
%nex=[14 15 16];
%nk=4;
[bigVAR, yy0, yy2] = vecm0(xx, xvec, beta, opts_vecm_.nex, opts_vecm_.nk);
[bigVAR1, yy10, yy12] = vecm2(xx, xvec, beta, opts_vecm_.nex, opts_vecm_.nk, 1);
save([fname_,'_vecm'],'yy0','yy2','bigVAR','yy10','yy12','bigVAR1')
% yy2 out of sample
% yy0 within sample (plus projections)
xxx=[xx xvec];
ifig=0;disp(' ')
%disp(' ')
%disp(' VECM RMSE')
for j=1:length(varargin),
rmseVECM(j)=sqrt(mean((xxx(istart:end,j)-squeeze(yy0(istart:size(xx,1),j,1))).^2));
rmseVECM1(j)=sqrt(mean((xxx(istart:end,j)-squeeze(yy10(istart:size(xx,1),j,1))).^2));
if gflag(j)
ttit=['diff(',varargin{j},')'];
else
ttit=varargin{j};
end
%disp([ttit, sprintf('%15.5g',[rmseVECM(j)'])])
if mod(j,9)==1,
figure('Name',['VECM 1-step ahead prediction']),
ifig=ifig+1;
end
subplot(3,3,j-9*(ifig-1))
plot(T(fobs+istart-1:fobs+nobs-1),[xxx(istart:end,j) squeeze(yy0(istart:size(xx,1),j,1)) squeeze(yy10(istart:size(xx,1),j,1))]) % example within sample plot
title(ttit,'interpreter','none')
if mod(j,9)==0 | j==length(varargin),
saveas(gcf,[fname_,'_VECM_Fit',int2str(ifig)])
if options_.nograph
close(gcf)
end
end
end
nk=opts_vecm_.nk;
ifig=0;
%disp(' ')
%disp(' ')
%disp(' VECM k-step RMSE')
for j=1:length(varargin),
rmseKVECM(j)=sqrt(mean((xxx(nk+1:end,j)-squeeze(yy0(nk+1:size(xx,1),j,nk))).^2));
rmseKVECM1(j)=sqrt(mean((xxx(nk+1:end,j)-squeeze(yy10(nk+1:size(xx,1),j,nk))).^2));
if gflag(j)
ttit=['diff(',varargin{j},')'];
else
ttit=varargin{j};
end
%disp([ttit, sprintf('%15.5g',[rmseKVECM(j)'])])
if mod(j,9)==1,
figure('Name',['VECM ',int2str(nk),'-step ahead prediction']),
ifig=ifig+1;
end
subplot(3,3,j-9*(ifig-1))
plot(T(fobs+nk:fobs+nobs-1),[xxx(nk+1:end,j) squeeze(yy0(nk+1:size(xx,1),j,nk)) squeeze(yy10(nk+1:size(xx,1),j,nk))]) % example within sample plot
title(ttit,'interpreter','none')
if mod(j,9)==0 | j==length(varargin),
saveas(gcf,[fname_,'_VECM_',num2str(nk),'-step',int2str(ifig)])
if options_.nograph
close(gcf)
end
end
end
% out of sample
% figure, plot([xxx(:,j)]),
% hold on, plot([size(xx,1):size(xx,1)+nk], [xxx(end,j); yy2(:,j)],'w')
% hold on, plot([size(xx,1):size(xx,1)+nk], [xxx(end,j); yy12(:,j)],'c')
% legend('data','projection')