1function [ConditionalVarianceDecomposition, ConditionalVarianceDecomposition_ME]= conditional_variance_decomposition(StateSpaceModel, Steps, SubsetOfVariables,sigma_e_is_diagonal)
2% This function computes the conditional variance decomposition of a given state space model
3% for a subset of endogenous variables.
4%
5% INPUTS
6%   StateSpaceModel     [structure]   Specification of the state space model.
7%   Steps               [integer]     1*h vector of dates.
8%   SubsetOfVariables   [integer]     1*q vector of indices (declaration order).
9%
10% OUTPUTS
11%   ConditionalVarianceDecomposition  [double] [n h p] array, where
12%                                                    n is equal to length(SubsetOfVariables)
13%                                                    h is the number of Steps
14%                                                    p is the number of state innovations and
15%   ConditionalVarianceDecomposition_ME  [double] [m h p] array, where
16%                                                    m is equal to length(intersect(SubsetOfVariables,varobs))
17%                                                    h is the number of Steps
18%                                                    p is the number of state innovations and
19
20% Copyright (C) 2010-2017 Dynare Team
21%
22% This file is part of Dynare.
23%
24% Dynare is free software: you can redistribute it and/or modify
25% it under the terms of the GNU General Public License as published by
26% the Free Software Foundation, either version 3 of the License, or
27% (at your option) any later version.
28%
29% Dynare is distributed in the hope that it will be useful,
30% but WITHOUT ANY WARRANTY; without even the implied warranty of
31% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
32% GNU General Public License for more details.
33%
34% You should have received a copy of the GNU General Public License
35% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
36
37if any(Steps <= 0)
38    error(['Conditional variance decomposition: All periods must be strictly ' ...
39           'positive'])
40end
41
42number_of_state_innovations = ...
43    StateSpaceModel.number_of_state_innovations;
44transition_matrix = StateSpaceModel.transition_matrix;
45number_of_state_equations = ...
46    StateSpaceModel.number_of_state_equations;
47order_var = StateSpaceModel.order_var;
48nSteps = length(Steps);
49
50ConditionalVariance = zeros(number_of_state_equations,nSteps,number_of_state_innovations);
51
52if StateSpaceModel.sigma_e_is_diagonal
53    B = StateSpaceModel.impulse_matrix.* ...
54        repmat(sqrt(diag(StateSpaceModel.state_innovations_covariance_matrix)'),...
55               number_of_state_equations,1);
56else
57    B = StateSpaceModel.impulse_matrix*chol(StateSpaceModel.state_innovations_covariance_matrix)';
58end
59
60for i=1:number_of_state_innovations
61    BB = B(:,i)*B(:,i)';
62    V = zeros(number_of_state_equations,number_of_state_equations);
63    m = 1;
64    for h = 1:max(Steps)
65        V = transition_matrix*V*transition_matrix'+BB;
66        if h == Steps(m)
67            ConditionalVariance(order_var,m,i) = diag(V);
68            m = m+1;
69        end
70    end
71end
72
73ConditionalVariance = ConditionalVariance(SubsetOfVariables,:,:);
74
75NumberOfVariables = length(SubsetOfVariables);
76SumOfVariances = zeros(NumberOfVariables,nSteps);
77for h = 1:length(Steps)
78    SumOfVariances(:,h) = sum(ConditionalVariance(:,h,:),3);
79end
80
81ConditionalVarianceDecomposition = zeros(NumberOfVariables,length(Steps),number_of_state_innovations);
82for i=1:number_of_state_innovations
83    for h = 1:length(Steps)
84        ConditionalVarianceDecomposition(:,h,i) = squeeze(ConditionalVariance(:,h,i))./SumOfVariances(:,h);
85    end
86end
87
88% get intersection of requested variables and observed variables with
89% Measurement error
90
91if ~all(diag(StateSpaceModel.measurement_error)==0)
92    if (isoctave && octave_ver_less_than('6')) || (~isoctave && matlab_ver_less_than('8.1'))
93        [observable_pos,index_subset,index_observables]=intersect_stable(SubsetOfVariables,StateSpaceModel.observable_pos);
94    else
95        [observable_pos,index_subset,index_observables]=intersect(SubsetOfVariables,StateSpaceModel.observable_pos,'stable');
96    end
97    ME_Variance=diag(StateSpaceModel.measurement_error);
98
99    ConditionalVarianceDecomposition_ME = zeros(length(observable_pos),length(Steps),number_of_state_innovations+1);
100    for i=1:number_of_state_innovations
101        for h = 1:length(Steps)
102            ConditionalVarianceDecomposition_ME(:,h,i) = squeeze(ConditionalVariance(index_subset,h,i))./(SumOfVariances(index_subset,h)+ME_Variance(index_observables));
103        end
104    end
105    ConditionalVarianceDecomposition_ME(:,:,number_of_state_innovations+1)=1-sum(ConditionalVarianceDecomposition_ME(:,:,1:number_of_state_innovations),3);
106else
107    ConditionalVarianceDecomposition_ME=[];
108end
109