From 1fc8bbd6d3dc3b280f28df247ebd45344b21a057 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Sun, 28 Nov 2010 11:03:40 +0100 Subject: [PATCH] bugs correction in computation of posterior moments for -conditional variance decomposition -hpdsup -moments with no variance in their posterior distribution modification of computation of conditional variance decomposition --- matlab/compute_moments_varendo.m | 2 +- matlab/conditional_variance_decomposition.m | 28 ++++++++------ ...ional_variance_decomposition_mc_analysis.m | 37 +++++++++---------- ...splay_conditional_variance_decomposition.m | 8 ++-- ...tical_conditional_variance_decomposition.m | 11 +++--- matlab/posterior_moments.m | 14 ++++--- 6 files changed, 55 insertions(+), 45 deletions(-) diff --git a/matlab/compute_moments_varendo.m b/matlab/compute_moments_varendo.m index a9ee76959..8db1c516e 100644 --- a/matlab/compute_moments_varendo.m +++ b/matlab/compute_moments_varendo.m @@ -113,7 +113,7 @@ if M_.exo_nbr > 1 if posterior for i=1:NumberOfEndogenousVariables for j=1:NumberOfExogenousVariables - oo_ = posterior_analysis('conditional decomposition',var_list_(i,:),M_.exo_names(j,:),Steps,options_,M_,oo_); + oo_ = posterior_analysis('conditional decomposition',i,M_.exo_names(j,:),Steps,options_,M_,oo_); end end else diff --git a/matlab/conditional_variance_decomposition.m b/matlab/conditional_variance_decomposition.m index bf8eef091..b9893ad4e 100644 --- a/matlab/conditional_variance_decomposition.m +++ b/matlab/conditional_variance_decomposition.m @@ -1,4 +1,4 @@ -function PackedConditionalVarianceDecomposition = conditional_variance_decomposition(StateSpaceModel, Steps, SubsetOfVariables,sigma_e_is_diagonal) +function ConditionalVarianceDecomposition = conditional_variance_decomposition(StateSpaceModel, Steps, SubsetOfVariables,sigma_e_is_diagonal) % This function computes the conditional variance decomposition of a given state space model % for a subset of endogenous variables. % @@ -8,9 +8,10 @@ function PackedConditionalVarianceDecomposition = conditional_variance_decomposi % SubsetOfVariables [integer] 1*q vector of indices. % % OUTPUTS -% PackedConditionalVarianceDecomposition [double] n(n+1)/2*p matrix, where p is the number of state innovations and -% n is equal to length(SubsetOfVariables). -% +% ConditionalVarianceDecomposition [double] [n h p] array, where +% n is equal to length(SubsetOfVariables) +% h is the number of Steps +% p is the number of state innovations and % SPECIAL REQUIREMENTS % % [1] In this version, absence of measurement errors is assumed... @@ -37,11 +38,10 @@ number_of_state_innovations = ... transition_matrix = StateSpaceModel.transition_matrix; number_of_state_equations = ... StateSpaceModel.number_of_state_equations; +order_var = StateSpaceModel.order_var; nSteps = length(Steps); -ConditionalVariance = zeros(number_of_state_equations,number_of_state_equations); -ConditionalVariance = repmat(ConditionalVariance,[1 1 nSteps ... - number_of_state_innovations]); +ConditionalVariance = zeros(number_of_state_equations,nSteps,number_of_state_innovations); if StateSpaceModel.sigma_e_is_diagonal B = StateSpaceModel.impulse_matrix.* ... @@ -58,17 +58,23 @@ for i=1:number_of_state_innovations for h = 1:max(Steps) V = transition_matrix*V*transition_matrix'+BB; if h == Steps(m) - ConditionalVariance(:,:,m,i) = V; + ConditionalVariance(order_var,m,i) = diag(V); m = m+1; end end end -ConditionalVariance = ConditionalVariance(SubsetOfVariables,SubsetOfVariables,:,:); +ConditionalVariance = ConditionalVariance(SubsetOfVariables,:,:); + NumberOfVariables = length(SubsetOfVariables); -PackedConditionalVarianceDecomposition = zeros(NumberOfVariables*(NumberOfVariables+1)/2,length(Steps),StateSpaceModel.number_of_state_innovations); +SumOfVariances = zeros(NumberOfVariables,nSteps); +for h = 1:length(Steps) + SumOfVariances(:,h) = sum(ConditionalVariance(:,h,:),3); +end + +ConditionalVarianceDecomposition = zeros(NumberOfVariables,length(Steps),number_of_state_innovations); for i=1:number_of_state_innovations for h = 1:length(Steps) - PackedConditionalVarianceDecomposition(:,h,i) = dyn_vech(ConditionalVariance(:,:,h,i)); + ConditionalVarianceDecomposition(:,h,i) = squeeze(ConditionalVariance(:,h,i))./SumOfVariances(:,h); end end \ No newline at end of file diff --git a/matlab/conditional_variance_decomposition_mc_analysis.m b/matlab/conditional_variance_decomposition_mc_analysis.m index 0968535f9..f96a67601 100644 --- a/matlab/conditional_variance_decomposition_mc_analysis.m +++ b/matlab/conditional_variance_decomposition_mc_analysis.m @@ -1,4 +1,5 @@ -function oo_ = conditional_variance_decomposition_mc_analysis(NumberOfSimulations, type, dname, fname, Steps, exonames, exo, vartan, var, mh_conf_sig, oo_) +function oo_ = ... + conditional_variance_decomposition_mc_analysis(NumberOfSimulations, type, dname, fname, Steps, exonames, exo, var_list, endogenous_variable_index, mh_conf_sig, oo_) % This function analyses the (posterior or prior) distribution of the % endogenous conditional variance decomposition. @@ -27,19 +28,19 @@ else PATH = [dname '/prior/moments/']; end -indx = check_name(vartan,var); -if isempty(indx) - disp([ type '_analysis:: ' var ' is not a stationary endogenous variable!']) - return -end -endogenous_variable_index = sum(1:indx); +% $$$ indx = check_name(vartan,var); +% $$$ if isempty(indx) +% $$$ disp([ type '_analysis:: ' var ' is not a stationary endogenous variable!']) +% $$$ return +% $$$ end +% $$$ endogenous_variable_index = sum(1:indx); exogenous_variable_index = check_name(exonames,exo); if isempty(exogenous_variable_index) disp([ type '_analysis:: ' exo ' is not a declared exogenous variable!']) return end -name = [ var '.' exo ]; +name = [ var_list(endogenous_variable_index,:) '.' exo ]; if isfield(oo_, [ TYPE 'TheoreticalMoments' ]) eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments;']) if isfield(temporary_structure,'dsge') @@ -75,17 +76,15 @@ p_density = NaN(2^9,2,length(Steps)); p_hpdinf = NaN(1,length(Steps)); p_hpdsup = NaN(1,length(Steps)); for i=1:length(Steps) - if ~isconst(tmp(:,i)) - [pp_mean, pp_median, pp_var, hpd_interval, pp_deciles, pp_density] = ... - posterior_moments(tmp(:,i),1,mh_conf_sig); - p_mean(2,i) = pp_mean; - p_median(i) = pp_median; - p_variance(i) = pp_var; - p_deciles(:,i) = pp_deciles; - p_hpdinf(i) = hpd_interval(1); - p_hpdinf(i) = hpd_interval(2); - p_density(:,:,i) = pp_density; - end + [pp_mean, pp_median, pp_var, hpd_interval, pp_deciles, pp_density] = ... + posterior_moments(tmp(:,i),1,mh_conf_sig); + p_mean(2,i) = pp_mean; + p_median(i) = pp_median; + p_variance(i) = pp_var; + p_deciles(:,i) = pp_deciles; + p_hpdinf(i) = hpd_interval(1); + p_hpdsup(i) = hpd_interval(2); + p_density(:,:,i) = pp_density; end eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.mean.' name ' = p_mean;']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.median.' name ' = p_median;']); diff --git a/matlab/display_conditional_variance_decomposition.m b/matlab/display_conditional_variance_decomposition.m index c4ad9b752..a75b95820 100644 --- a/matlab/display_conditional_variance_decomposition.m +++ b/matlab/display_conditional_variance_decomposition.m @@ -44,8 +44,9 @@ ic = dr.nstatic+(1:dr.npred)'; [StateSpaceModel.transition_matrix,StateSpaceModel.impulse_matrix] = kalman_transition_matrix(dr,iv,ic,[],exo_nbr); StateSpaceModel.state_innovations_covariance_matrix = M_.Sigma_e; +StateSpaceModel.order_var = dr.order_var; -conditional_decomposition_array = conditional_variance_decomposition(StateSpaceModel,Steps,dr.inv_order_var(SubsetOfVariables )); +conditional_decomposition_array = conditional_variance_decomposition(StateSpaceModel,Steps,SubsetOfVariables ); if options_.noprint == 0 disp(' ') @@ -58,10 +59,9 @@ for i=1:length(Steps) disp(['Period ' int2str(Steps(i)) ':']) for j=1:exo_nbr - vardec_i(:,j) = dyn_diag_vech(conditional_decomposition_array(:, ... - i,j)); + vardec_i(:,j) = 100*conditional_decomposition_array(:, ... + i,j); end - vardec_i = 100*vardec_i./repmat(sum(vardec_i,2),1,exo_nbr); if options_.noprint == 0 headers = M_.exo_names; headers(M_.exo_names_orig_ord,:) = headers; diff --git a/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m b/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m index b35166dcb..f7543882f 100644 --- a/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m +++ b/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m @@ -67,14 +67,14 @@ nar = options_.ar; options_.ar = 0; NumberOfDrawsFiles = rows(DrawsFiles); -NumberOfSavedElementsPerSimulation = nvar*(nvar+1)/2*M_.exo_nbr*length(Steps); +NumberOfSavedElementsPerSimulation = nvar*M_.exo_nbr*length(Steps); MaXNumberOfConditionalDecompLines = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation/8); if SampleSize<=MaXNumberOfConditionalDecompLines - Conditional_decomposition_array = zeros(nvar*(nvar+1)/2,length(Steps),M_.exo_nbr,SampleSize); + Conditional_decomposition_array = zeros(nvar,length(Steps),M_.exo_nbr,SampleSize); NumberOfConditionalDecompFiles = 1; else - Conditional_decomposition_array = zeros(nvar*(nvar+1)/2,length(Steps),M_.exo_nbr,MaXNumberOfConditionalDecompLines); + Conditional_decomposition_array = zeros(nvar,length(Steps),M_.exo_nbr,MaXNumberOfConditionalDecompLines); NumberOfLinesInTheLastConditionalDecompFile = mod(SampleSize,MaXNumberOfConditionalDecompLines); NumberOfConditionalDecompFiles = ceil(SampleSize/MaXNumberOfConditionalDecompLines); end @@ -118,6 +118,7 @@ for file = 1:NumberOfDrawsFiles StateSpaceModel.number_of_state_equations = M_.endo_nbr+rows(aux); StateSpaceModel.number_of_state_innovations = M_.exo_nbr; StateSpaceModel.sigma_e_is_diagonal = M_.sigma_e_is_diagonal; + StateSpaceModel.order_var = dr.order_var; first_call = 0; clear('endo_nbr','nstatic','npred','k'); end @@ -136,10 +137,10 @@ for file = 1:NumberOfDrawsFiles 'Conditional_decomposition_array'); end if (ConditionalDecompFileNumber==NumberOfConditionalDecompFiles-1)% Prepare last round. - Conditional_decomposition_array = zeros(nvar*(nvar+1)/2, length(Steps),M_.exo_nbr,NumberOfLinesInTheLastConditionalDecompFile) ; + Conditional_decomposition_array = zeros(nvar, length(Steps),M_.exo_nbr,NumberOfLinesInTheLastConditionalDecompFile) ; NumberOfConditionalDecompLines = NumberOfLinesInTheLastConditionalDecompFile; elseif ConditionalDecompFileNumber 0 + bandwidth = 0; % Rule of thumb optimal bandwidth parameter. + kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton. + optimal_bandwidth = mh_optimal_bandwidth(xx,number_of_draws,bandwidth,kernel_function); + [density(:,1),density(:,2)] = kernel_density_estimate(xx,number_of_grid_points,... + number_of_draws,optimal_bandwidth,kernel_function); + else + density = NaN(number_of_grid_points,2); + end end \ No newline at end of file