dsge_posterior_kernel.m 11.4 KB
Newer Older
1
2
function [fval,cost_flag,ys,trend_coeff,info] = dsge_posterior_kernel(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
% function [fval,cost_flag,ys,trend_coeff,info] = dsge_posterior_kernel(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
% Evaluates the posterior kernel of a dsge model. 
% 
% INPUTS 
%   xparam1                        [double]   vector of model parameters.
%   gend                           [integer]  scalar specifying the number of observations.
%   data                           [double]   matrix of data
%   data_index                     [cell]     cell of column vectors
%   number_of_observations         [integer]
%   no_more_missing_observations   [integer] 
% OUTPUTS 
%   fval        :     value of the posterior kernel at xparam1.
%   cost_flag   :     zero if the function returns a penalty, one otherwise.
%   ys          :     steady state of original endogenous variables
%   trend_coeff :
%   info        :     vector of informations about the penalty:
%                     41: one (many) parameter(s) do(es) not satisfied the lower bound
%                     42: one (many) parameter(s) do(es) not satisfied the upper bound
%               
% SPECIAL REQUIREMENTS
%

24
% Copyright (C) 2004-2011 Dynare Team
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.

41
global bayestopt_ estim_params_ options_ trend_coeff_ M_ oo_
42
43
44
45
46
fval            = [];
ys              = [];
trend_coeff     = [];
cost_flag       = 1;
nobs            = size(options_.varobs,1);
47
48
49
%------------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
50
if ~isequal(options_.mode_compute,1) && any(xparam1 < bayestopt_.lb)
51
52
53
54
55
56
    k = find(xparam1 < bayestopt_.lb);
    fval = bayestopt_.penalty+sum((bayestopt_.lb(k)-xparam1(k)).^2);
    cost_flag = 0;
    info = 41;
    return;
end
57
if ~isequal(options_.mode_compute,1) && any(xparam1 > bayestopt_.ub)
58
59
60
61
62
63
64
65
66
    k = find(xparam1 > bayestopt_.ub);
    fval = bayestopt_.penalty+sum((xparam1(k)-bayestopt_.ub(k)).^2);
    cost_flag = 0;
    info = 42;
    return;
end
Q = M_.Sigma_e;
H = M_.H;
for i=1:estim_params_.nvx
67
68
    k =estim_params_.var_exo(i,1);
    Q(k,k) = xparam1(i)*xparam1(i);
69
70
71
end
offset = estim_params_.nvx;
if estim_params_.nvn
72
    for i=1:estim_params_.nvn
73
74
        k = estim_params_.var_endo(i,1);
        H(k,k) = xparam1(i+offset)*xparam1(i+offset);
75
76
    end
    offset = offset+estim_params_.nvn;
77
end     
78
if estim_params_.ncx
79
    for i=1:estim_params_.ncx
80
81
82
83
        k1 =estim_params_.corrx(i,1);
        k2 =estim_params_.corrx(i,2);
        Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2));
        Q(k2,k1) = Q(k1,k2);
84
85
    end
    [CholQ,testQ] = chol(Q);
86
    if testQ    %% The variance-covariance matrix of the structural innovations is not definite positive.
87
88
89
90
91
92
93
94
95
        %% We have to compute the eigenvalues of this matrix in order to build the penalty.
        a = diag(eig(Q));
        k = find(a < 0);
        if k > 0
            fval = bayestopt_.penalty+sum(-a(k));
            cost_flag = 0;
            info = 43;
            return
        end
96
97
    end
    offset = offset+estim_params_.ncx;
98
99
end
if estim_params_.ncn 
100
    for i=1:estim_params_.ncn
101
102
103
104
        k1 = options_.lgyidx2varobs(estim_params_.corrn(i,1));
        k2 = options_.lgyidx2varobs(estim_params_.corrn(i,2));
        H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
        H(k2,k1) = H(k1,k2);
105
106
107
    end
    [CholH,testH] = chol(H);
    if testH
108
109
110
111
112
113
114
115
        a = diag(eig(H));
        k = find(a < 0);
        if k > 0
            fval = bayestopt_.penalty+sum(-a(k));
            cost_flag = 0;
            info = 44;
            return
        end
116
117
    end
    offset = offset+estim_params_.ncn;
118
119
120
121
122
123
124
125
126
end
if estim_params_.np > 0
    M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end);
end
M_.Sigma_e = Q;
M_.H = H;
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
127
[T,R,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
128
if info(1) == 1 || info(1) == 2 || info(1) == 5
129
130
131
    fval = bayestopt_.penalty+1;
    cost_flag = 0;
    return
132
elseif info(1) == 3 || info(1) == 4 || info(1) == 20
133
134
135
    fval = bayestopt_.penalty+info(2);%^2; % penalty power raised in DR1.m and resol already. GP July'08
    cost_flag = 0;
    return
136
137
138
end
bayestopt_.mf = bayestopt_.mf1;
if ~options_.noconstant
139
    if options_.loglinear == 1
140
        constant = log(SteadyState(bayestopt_.mfys));
141
    else
142
        constant = SteadyState(bayestopt_.mfys);
143
    end
144
else
145
    constant = zeros(nobs,1);
146
147
end
if bayestopt_.with_trend == 1
148
149
150
    trend_coeff = zeros(nobs,1);
    t = options_.trend_coeffs;
    for i=1:length(t)
151
152
153
        if ~isempty(t{i})
            trend_coeff(i) = evalin('base',t{i});
        end
154
155
    end
    trend = repmat(constant,1,gend)+trend_coeff*[1:gend];
156
else
157
    trend = repmat(constant,1,gend);
158
159
160
161
162
163
164
165
166
end
start = options_.presample+1;
np    = size(T,1);
mf    = bayestopt_.mf;
no_missing_data_flag = (number_of_observations==gend*nobs);
%------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
%------------------------------------------------------------------------------
kalman_algo = options_.kalman_algo;
167
if options_.lik_init == 1               % Kalman filter
168
169
170
171
    if kalman_algo ~= 2
        kalman_algo = 1;
    end
    Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);
172
173
    Pinf        = [];
elseif options_.lik_init == 2   % Old Diffuse Kalman filter
174
175
176
177
178
    if kalman_algo ~= 2
        kalman_algo = 1;
    end
    Pstar = 10*eye(np);
    Pinf = [];
179
elseif options_.lik_init == 3   % Diffuse Kalman filter
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
    if kalman_algo ~= 4
        kalman_algo = 3;
    end
    [QT,ST] = schur(T);
    e1 = abs(ordeig(ST)) > 2-options_.qz_criterium;
    [QT,ST] = ordschur(QT,ST,e1);
    k = find(abs(ordeig(ST)) > 2-options_.qz_criterium);
    nk = length(k);
    nk1 = nk+1;
    Pinf = zeros(np,np);
    Pinf(1:nk,1:nk) = eye(nk);
    Pstar = zeros(np,np);
    B = QT'*R*Q*R'*QT;
    for i=np:-1:nk+2
        if ST(i,i-1) == 0
            if i == np
                c = zeros(np-nk,1);
            else
                c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
                    ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
            end
            q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i);
            Pstar(nk1:i,i) = q\(B(nk1:i,i)+c);
            Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
        else
            if i == np
                c = zeros(np-nk,1);
                c1 = zeros(np-nk,1);
            else
                c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
                    ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i)+...
                    ST(i,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1);
                c1 = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i-1,i+1:end)')+...
                     ST(i-1,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1)+...
                     ST(i-1,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
            end
            q = [eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i) -ST(nk1:i,nk1:i)*ST(i,i-1);...
                 -ST(nk1:i,nk1:i)*ST(i-1,i) eye(i-nk)-ST(nk1:i,nk1:i)*ST(i-1,i-1)];
            z =  q\[B(nk1:i,i)+c;B(nk1:i,i-1)+c1];
            Pstar(nk1:i,i) = z(1:(i-nk));
            Pstar(nk1:i,i-1) = z(i-nk+1:end);
            Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
            Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)';
            i = i - 1;
        end
    end
    if i == nk+2
        c = ST(nk+1,:)*(Pstar(:,nk+2:end)*ST(nk1,nk+2:end)')+ST(nk1,nk1)*ST(nk1,nk+2:end)*Pstar(nk+2:end,nk1);
        Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1));
    end
    Z = QT(mf,:);
    R1 = QT'*R;
    [QQ,RR,EE] = qr(Z*ST(:,1:nk),0);
    k = find(abs(diag(RR)) < 1e-8);
    if length(k) > 0
        k1 = EE(:,k);
        dd =ones(nk,1);
        dd(k1) = zeros(length(k1),1);
        Pinf(1:nk,1:nk) = diag(dd);
    end
end
if kalman_algo == 2
    no_correlation_flag = 1;
    if length(H)==1
        H = zeros(nobs,1);
    else
        if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
            H = diag(H);
        else
            no_correlation_flag = 0;
        end
    end
end
kalman_tol = options_.kalman_tol;
riccati_tol = options_.riccati_tol;
mf = bayestopt_.mf1;
Y   = data-trend;
%------------------------------------------------------------------------------
% 4. Likelihood evaluation
%------------------------------------------------------------------------------
if (kalman_algo==1)% Multivariate Kalman Filter
    if no_missing_data_flag
        LIK = kalman_filter(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol); 
    else
        LIK = ...
            missing_observations_kalman_filter(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol, ...
                                               data_index,number_of_observations,no_more_missing_observations);
    end
    if isinf(LIK)
        kalman_algo = 2;
    end
end
if (kalman_algo==2)% Univariate Kalman Filter
    if no_correlation_flag
        LIK = univariate_kalman_filter(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol,data_index,number_of_observations,no_more_missing_observations);
    else
        LIK = univariate_kalman_filter_corr(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol,data_index,number_of_observations,no_more_missing_observations);
    end
end
if (kalman_algo==3)% Multivariate Diffuse Kalman Filter
    if no_missing_data_flag
        LIK = diffuse_kalman_filter(ST,R1,Q,H,Pinf,Pstar,Y,start,Z,kalman_tol,riccati_tol);
    else
        LIK = missing_observations_diffuse_kalman_filter(ST,R1,Q,H,Pinf,Pstar,Y,start,Z,kalman_tol,riccati_tol,...
                                                         data_index,number_of_observations,no_more_missing_observations);
    end
    if isinf(LIK)
        kalman_algo = 4;
    end
end
if (kalman_algo==4)% Univariate Diffuse Kalman Filter
    if no_correlation_flag
        LIK = univariate_diffuse_kalman_filter(ST,R1,Q,H,Pinf,Pstar,Y,start,Z,kalman_tol,riccati_tol,...
                                               data_index,number_of_observations,no_more_missing_observations);
    else
        LIK = univariate_diffuse_kalman_filter_corr(ST,R1,Q,H,Pinf,Pstar,Y,start,Z,kalman_tol,riccati_tol,...
                                                    data_index,number_of_observations,no_more_missing_observations);
    end
end
if imag(LIK) ~= 0
    likelihood = bayestopt_.penalty;
else
    likelihood = LIK;
end
% ------------------------------------------------------------------------------
% Adds prior if necessary
% ------------------------------------------------------------------------------
lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
fval    = (likelihood-lnprior);
309
options_.kalman_algo = kalman_algo;