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

@@ -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

matlab/conditional_variance_decomposition.m

-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

matlab/conditional_variance_decomposition_mc_analysis.m

-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;']);

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;

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<NumberOfConditionalDecompFiles-1
-                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);
             else
                 clear('Conditional_decomposition_array');
             end

matlab/posterior_moments.m

@@ -70,9 +70,13 @@ post_deciles = xx([round(0.1*number_of_draws) ...
 density = [];
 if info
     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.  
-    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);
+    if post_var > 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