DsgeVarLikelihood.m 8.77 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
function [fval,cost_flag,info,PHI,SIGMAu,iXX,prior] = DsgeVarLikelihood(xparam1,gend)
% Evaluates the posterior kernel of the bvar-dsge model. 
% 
% INPUTS 
%   o xparam1       [double]     Vector of model's parameters.
%   o gend          [integer]    Number of observations (without conditionning observations for the lags).
%  
% OUTPUTS 
%   o fval          [double]     Value of the posterior kernel at xparam1.
%   o cost_flag     [integer]    Zero if the function returns a penalty, one otherwise.
%   o info          [integer]    Vector of informations about the penalty.
%   o PHI           [double]     Stacked BVAR-DSGE autoregressive matrices (at the mode associated to xparam1).
%   o SIGMAu        [double]     Covariance matrix of the BVAR-DSGE (at the mode associated to xparam1).
%   o iXX           [double]     inv(X'X).
%   o prior         [double]     a matlab structure describing the dsge-var prior.  
%
% SPECIAL REQUIREMENTS
%   None.

20
% Copyright (C) 2006-2011 Dynare Team
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
%
% 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/>.

37
global bayestopt_ estim_params_ M_ options_ oo_
38
39
40
41
42
43
44
45
46
47

nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np  = estim_params_.np;
nx = nvx+nvn+ncx+ncn+np;
ns = nvx+nvn+ncx+ncn;

NumberOfObservedVariables = size(options_.varobs,1);
48
NumberOfLags = options_.dsge_varlag;
49
50
51
52
53
54
55
56
57
58
59
60
61
NumberOfParameters = NumberOfObservedVariables*NumberOfLags ;
if ~options_.noconstant
    NumberOfParameters = NumberOfParameters + 1;
end

mYY = evalin('base', 'mYY');
mYX = evalin('base', 'mYX');
mXY = evalin('base', 'mXY');
mXX = evalin('base', 'mXX');

fval = [];
cost_flag = 1;

62
if ~isequal(options_.mode_compute,1) && any(xparam1 < bayestopt_.lb)
63
64
65
66
67
68
69
    k = find(xparam1 < bayestopt_.lb);
    fval = bayestopt_.penalty+sum((bayestopt_.lb(k)-xparam1(k)).^2);
    cost_flag = 0;
    info = 41;
    return;
end

70
if ~isequal(options_.mode_compute,11) && any(xparam1 > bayestopt_.ub)
71
72
73
74
75
76
77
78
79
80
81
82
83
84
    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;
for i=1:estim_params_.nvx
    k = estim_params_.var_exo(i,1);
    Q(k,k) = xparam1(i)*xparam1(i);
end
offset = estim_params_.nvx;
if estim_params_.nvn
85
    disp('DsgeVarLikelihood :: Measurement errors are not implemented!')
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
    return
end 
if estim_params_.ncx
    disp('DsgeVarLikelihood :: Correlated structural innovations are not implemented!')
    return
end

M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end);
M_.Sigma_e = Q;

%% Weight of the dsge prior:
dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names));
% Is the DSGE prior proper?
if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/gend;
    fval = bayestopt_.penalty+abs(gend*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables));
    cost_flag = 0;
    info = 51;
    return;
end

%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
109
[T,R,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
110

111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
if info(1) == 1 || info(1) == 2 || info(1) == 5
    fval = bayestopt_.penalty+1;
    cost_flag = 0;
    return
elseif info(1) == 3 || info(1) == 4 || info(1) == 19 || info(1) == 20 || info(1) == 21
    fval = bayestopt_.penalty+info(2);
    cost_flag = 0;
    return
end

if ~options_.noconstant
    if options_.loglinear 
        constant = transpose(log(SteadyState(bayestopt_.mfys)));
    else
        constant = transpose(SteadyState(bayestopt_.mfys));
126
    end 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
else
    constant = zeros(1,NumberOfObservedVariables);
end
if bayestopt_.with_trend == 1
    disp('DsgeVarLikelihood :: Linear trend is not yet implemented!')
    return
end

%------------------------------------------------------------------------------
% 3. theoretical moments (second order)
%------------------------------------------------------------------------------
tmp0 = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);% I compute the variance-covariance matrix
mf  = bayestopt_.mf1;          % of the restricted state vector.

% Get the non centered second order moments
TheoreticalAutoCovarianceOfTheObservedVariables = ...
    zeros(NumberOfObservedVariables,NumberOfObservedVariables,NumberOfLags+1);
TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) = tmp0(mf,mf)+constant'*constant;
for lag = 1:NumberOfLags
146
147
148
    tmp0 = T*tmp0;
    TheoreticalAutoCovarianceOfTheObservedVariables(:,:,lag+1) = tmp0(mf,mf) ...
        + constant'*constant;
149
150
151
152
end
% Build the theoretical "covariance" between Y and X
GYX = zeros(NumberOfObservedVariables,NumberOfParameters);
for i=1:NumberOfLags
153
154
    GYX(:,(i-1)*NumberOfObservedVariables+1:i*NumberOfObservedVariables) = ...
        TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1);
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
end
if ~options_.noconstant
    GYX(:,end) = constant';
end
% Build the theoretical "covariance" between X and X
GXX = kron(eye(NumberOfLags), ...
           TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1));
for i = 1:NumberOfLags-1
    tmp1 = diag(ones(NumberOfLags-i,1),i); 
    tmp2 = diag(ones(NumberOfLags-i,1),-i);
    GXX = GXX + kron(tmp1,TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1));
    GXX = GXX + kron(tmp2,TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1)');
end

if ~options_.noconstant
    % Add one row and one column to GXX
    GXX = [GXX , kron(ones(NumberOfLags,1),constant') ; [  kron(ones(1,NumberOfLags),constant) , 1] ];
end

GYY = TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1);

assignin('base','GYY',GYY);
assignin('base','GXX',GXX);
assignin('base','GYX',GYX);

if ~isinf(dsge_prior_weight)
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
    tmp0 = dsge_prior_weight*gend*TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) + mYY ;
    tmp1 = dsge_prior_weight*gend*GYX + mYX;
    tmp2 = inv(dsge_prior_weight*gend*GXX+mXX);
    SIGMAu = tmp0 - tmp1*tmp2*tmp1'; clear('tmp0');
    if ~ispd(SIGMAu)
        v = diag(SIGMAu);
        k = find(v<0);
        fval = bayestopt_.penalty + sum(v(k).^2);
        info = 52;
        cost_flag = 0;
        return;
    end
    SIGMAu = SIGMAu / (gend*(1+dsge_prior_weight));
    PHI = tmp2*tmp1'; clear('tmp1');
    prodlng1 = sum(gammaln(.5*((1+dsge_prior_weight)*gend- ...
                               NumberOfObservedVariables*NumberOfLags ...
                               +1-(1:NumberOfObservedVariables)')));
    prodlng2 = sum(gammaln(.5*(dsge_prior_weight*gend- ...
                               NumberOfObservedVariables*NumberOfLags ...
                               +1-(1:NumberOfObservedVariables)')));  
    lik = .5*NumberOfObservedVariables*log(det(dsge_prior_weight*gend*GXX+mXX)) ...
          + .5*((dsge_prior_weight+1)*gend-NumberOfParameters)*log(det((dsge_prior_weight+1)*gend*SIGMAu)) ...
          - .5*NumberOfObservedVariables*log(det(dsge_prior_weight*gend*GXX)) ...
          - .5*(dsge_prior_weight*gend-NumberOfParameters)*log(det(dsge_prior_weight*gend*(GYY-GYX*inv(GXX)*GYX'))) ...
          + .5*NumberOfObservedVariables*gend*log(2*pi)  ...
          - .5*log(2)*NumberOfObservedVariables*((dsge_prior_weight+1)*gend-NumberOfParameters) ...
          + .5*log(2)*NumberOfObservedVariables*(dsge_prior_weight*gend-NumberOfParameters) ...
          - prodlng1 + prodlng2;
209
else
210
211
212
213
214
215
    iGXX = inv(GXX);
    SIGMAu = GYY - GYX*iGXX*transpose(GYX);
    PHI = iGXX*transpose(GYX);
    lik = gend * ( log(det(SIGMAu)) + NumberOfObservedVariables*log(2*pi) +  ...
                   trace(inv(SIGMAu)*(mYY - transpose(mYX*PHI) - mYX*PHI + transpose(PHI)*mXX*PHI)/gend));
    lik = .5*lik;% Minus likelihood
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
end      

lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
fval = (lik-lnprior);

if (nargout == 6)
    if isinf(dsge_prior_weight)
        iXX = iGXX;
    else
        iXX = tmp2;
    end
end

if (nargout==7)
    if isinf(dsge_prior_weight)
        iXX = iGXX;
    else
        iXX = tmp2;
    end
    iGXX = inv(GXX); 
    prior.SIGMAstar = GYY - GYX*iGXX*GYX';
    prior.PHIstar = iGXX*transpose(GYX);
    prior.ArtificialSampleSize = fix(dsge_prior_weight*gend);
    prior.DF = prior.ArtificialSampleSize - NumberOfParameters - NumberOfObservedVariables;
    prior.iGXX = iGXX;
241
end