Commit 4558db92 authored by rattoma's avatar rattoma
Browse files

Fixes to match new v4.1 changes.

git-svn-id: file:///home/sebastien/dynare/gsa_dyn@50 f1850c17-3b45-254b-b221-fcb05880fee1
parent 9378db0f
......@@ -53,6 +53,7 @@ rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet
options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1);
gend = options_.nobs;
n_varobs = size(options_.varobs,1);
rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
if options_.loglinear == 1 & ~options_.logdata
......@@ -64,6 +65,7 @@ if options_.prefilter == 1
else
data = transpose(rawdata);
end
[data_index,number_of_observations,no_more_missing_observations] = describe_missing_data(data,gend,n_varobs);
if ~isreal(rawdata)
error(['There are complex values in the data. Probably a wrong' ...
......@@ -89,7 +91,7 @@ load([OutDir,'\',namfile],'lpmat', 'lpmat0', 'istable')
x=[lpmat0(istable,:) lpmat(istable,:)];
clear lpmat lpmat0 istable %iunstable egg yys T
B = size(x,1);
[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff, aK] = DsgeSmoother(x(1,:)',gend,data);
[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff, aK] = DsgeSmoother(x(1,:)',gend,data,{},0);
n1=size(atT,1);
nfil=B/40;
......@@ -103,18 +105,19 @@ delete([OutDir,'\',namfile,'_*.mat'])
ib=0;
ifil=0;
opt_gsa=options_.opt_gsa;
options_.filter_step_ahead=1;
for b=1:B
ib=ib+1;
deep = x(b,:)';
set_all_parameters(deep);
dr = resol(oo_.steady_state,0);
%deep(1:offset) = xparam1(1:offset);
logpo2(b,1) = DsgeLikelihood(deep,gend,data);
logpo2(b,1) = DsgeLikelihood(deep,gend,data,data_index,number_of_observations,no_more_missing_observations);
if opt_gsa.lik_only==0,
[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff, aK] = DsgeSmoother(deep,gend,data);
[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff, aK] = DsgeSmoother(deep,gend,data,{},0);
stock_smooth(:,:,ib)=atT(1:M_.endo_nbr,:);
stock_filter(:,:,ib)=filtered_state_vector(1:M_.endo_nbr,:);
%stock_filter(:,:,ib)=aK(options_.filter_step_ahead,1:M_.endo_nbr,:);
% stock_filter(:,:,ib)=filtered_state_vector(1:M_.endo_nbr,:);
stock_filter(:,:,ib)=aK(1,1:M_.endo_nbr,:);
stock_ys(ib,:)=ys';
if ib==40,
ib=0;
......
......@@ -119,7 +119,7 @@ if ~loadSA,
else
%load([DirectoryName '/' M_.fname '_data.mat']);
[stock_gend, stock_data] = read_data;
filfilt = dir([DirectoryName '/' M_.fname '_filter*.mat']);
filfilt = dir([DirectoryName '/' M_.fname '_filter_step_ahead*.mat']);
filparam = dir([DirectoryName '/' M_.fname '_param*.mat']);
x=[];
logpo2=[];
......@@ -154,21 +154,23 @@ if ~loadSA,
jxj = strmatch(vj,lgy_(dr_.order_var,:),'exact');
js = strmatch(vj,lgy_,'exact');
if exist('xparam1','var')
if isfield(oo_,'FilteredVariables')
eval(['rmse_mode(i) = sqrt(mean((vobs(istart:end)-oo_.steady_state(js)-oo_.FilteredVariables.',vj,'(istart:end-1)).^2));'])
else
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1,stock_gend,stock_data);
y0 = ahat(jxj,:)' + ...
kron(ys_mode(js,:),ones(size(ahat,2),1));
% if isfield(oo_,'FilteredVariables')
% eval(['rmse_mode(i) = sqrt(mean((vobs(istart:end)-oo_.steady_state(js)-oo_.FilteredVariables.',vj,'(istart:end-1)).^2));'])
% else
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1,stock_gend,stock_data,{},0);
y0 = squeeze(aK(1,jxj,:)) + ...
kron(ys_mode(js,:),ones(size(aK,3),1));
% y0 = ahat(jxj,:)' + ...
% kron(ys_mode(js,:),ones(size(ahat,2),1));
rmse_mode(i) = sqrt(mean((vobs(istart:end)-y0(istart:end-1)).^2));
end
% end
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']);
load([DirectoryName '/' M_.fname '_filter_step_ahead',num2str(j),'.mat']);
nb = size(stock,4);
% y0(:,nbb+1:nbb+nb)=squeeze(stock(1,js,:,:)) + ...
% kron(sto_ys(nbb+1:nbb+nb,js)',ones(size(stock,3),1));
......@@ -200,9 +202,11 @@ if ~loadSA,
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));
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1_mean,stock_gend,stock_data,{},0);
y0 = squeeze(aK(1,jxj,:)) + ...
kron(ys_mean(js,:),ones(size(aK,3),1));
% y0 = ahat(jxj,:)' + ...
% kron(ys_mean(js,:),ones(size(ahat,2),1));
rmse_pmean(i) = sqrt(mean((vobs(istart:end)-y0(istart:end-1)).^2));
end
end
......@@ -554,7 +558,7 @@ for ix=1:ceil(length(nsnam)/5),
subplot(2,3,j-5*(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);
size(x,1),optimal_bandwidth,kernel_function);
h0 = plot(x1, f1,'k');
hold on,
np=find(SP(nsnam(j),:));
......@@ -563,7 +567,7 @@ for ix=1:ceil(length(nsnam)/5),
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);
nfilt, optimal_bandwidth,kernel_function);
h0 = plot(x1, f1);
set(h0,'color',a0(i,:))
end
......@@ -613,11 +617,11 @@ close all
% 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);
% nfilt, 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);
% nruns-nfilt,optimal_bandwidth,kernel_function);
% hold on, plot(x1, f1,'k','linewidth',2)
% ydum=get(gca,'ylim');
% %xdum=xparam1(nshock+np(i));
......
function pdraw = prior_draw_gsa(init,rdraw,cc)
function pdraw = prior_draw_gsa(init,rdraw)
% Draws from the prior distributions
% Adapted by M. Ratto from prior_draw (of DYNARE, copyright M. Juillard),
% for use with Sensitivity Toolbox for DYNARE
......@@ -7,9 +7,6 @@ function pdraw = prior_draw_gsa(init,rdraw,cc)
% INPUTS
% o init [integer] scalar equal to 1 (first call) or 0.
% o rdraw
% 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.
......@@ -37,65 +34,17 @@ function pdraw = prior_draw_gsa(init,rdraw,cc)
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
global M_ options_ estim_params_ bayestopt_
persistent fname npar bounds pshape pmean pstd a b p1 p2 p3 p4 condition
% global M_ options_ estim_params_ bayestopt_
global bayestopt_
persistent npar pshape p6 p7 p3 p4
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;
p6 = bayestopt_.p6;
p7 = bayestopt_.p7;
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;
npar = length(p6);
pdraw = zeros(npar,1);
return
end
......@@ -106,17 +55,17 @@ for i = 1:npar
case 5% Uniform prior.
pdraw(:,i) = rdraw(:,i)*(p4(i)-p3(i)) + p3(i);
case 3% Gaussian prior.
pdraw(:,i) = norm_inv(rdraw(:,i),pmean(i),pstd(i));
pdraw(:,i) = norm_inv(rdraw(:,i),p6(i),p7(i));
case 2% Gamma prior.
pdraw(:,i) = gamm_inv(rdraw(:,i),a(i),b(i))+p3(i);
pdraw(:,i) = gamm_inv(rdraw(:,i),p6(i),p7(i))+p3(i);
case 1% Beta distribution (TODO: generalized beta distribution)
pdraw(:,i) = beta_inv(rdraw(:,i),a(i),b(i))*(p4(i)-p3(i))+p3(i);
pdraw(:,i) = beta_inv(rdraw(:,i),p6(i),p7(i))*(p4(i)-p3(i))+p3(i);
case 4% INV-GAMMA1 distribution
% TO BE CHECKED
pdraw(:,i) = sqrt(1./gamm_inv(rdraw(:,i),p2(i)/2,2/p1(i)));
pdraw(:,i) = sqrt(1./gamm_inv(rdraw(:,i),p7(i)/2,2/p6(i)))+p3(i);
case 6% INV-GAMMA2 distribution
% TO BE CHECKED
pdraw(:,i) = 1./gamm_inv(rdraw(:,i),p2(i)/2,2/p1(i));
pdraw(:,i) = 1./gamm_inv(rdraw(:,i),p7(i)/2,2/p6(i))+p3(i);
otherwise
% Nothing to do here.
end
......
......@@ -123,7 +123,7 @@ if fload==0,
lpmat(:,j)=lpmat(randperm(Nsam),j);
end
end
else ilptau==0
else %ilptau==0
%[lpmat] = rand(Nsam,np);
for j=1:np,
lpmat(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
......@@ -131,16 +131,16 @@ if fload==0,
end
end
try
% try
dummy=prior_draw_gsa(1);
catch
if pprior,
if opt_gsa.prior_range==0;
error('Some unknown prior is specified or ML estimation,: use prior_range=1 option!!');
end
end
end
% catch
% if pprior,
% if opt_gsa.prior_range==0;
% error('Some unknown prior is specified or ML estimation,: use prior_range=1 option!!');
% end
% end
%
% end
if pprior,
for j=1:nshock,
if opt_gsa.morris~=1,
......@@ -151,11 +151,27 @@ if fload==0,
end
end
if opt_gsa.prior_range
if opt_gsa.identification,
deltx=min(0.001, 1/Nsam/2);
for j=1:np,
xdelt(:,:,j)=prior_draw_gsa(0,[lpmat0 lpmat]+deltx);
end
end
for j=1:np,
lpmat(:,j)=lpmat(:,j).*(bayestopt_.ub(j+nshock)-bayestopt_.lb(j+nshock))+bayestopt_.lb(j+nshock);
end
else
xx=prior_draw_gsa(0,[lpmat0 lpmat]);
if opt_gsa.identification,
deltx=min(0.001, 1/Nsam/2);
ldum=[lpmat0 lpmat];
ldum = prior_draw_gsa(0,ldum+deltx);
for j=1:nshock+np,
xdelt(:,:,j)=xx;
xdelt(:,j,j)=ldum(:,j);
end
clear ldum
end
lpmat0=xx(:,1:nshock);
lpmat=xx(:,nshock+1:end);
clear xx;
......
......@@ -19,7 +19,7 @@ function [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list)
end
end
[gamma_y,ivar] = th_autocovariances(dr,ivar);
[gamma_y,ivar] = th_autocovariances(dr,ivar,M_, options_);
m = dr.ys(ivar);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment