DsgeLikelihood.m 12.2 KB
Newer Older
1
function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
2
% function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
assia's avatar
assia committed
3 4 5
% Evaluates the posterior kernel of a dsge model. 
% 
% INPUTS 
6 7 8 9 10 11
%   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] 
assia's avatar
assia committed
12
% OUTPUTS 
13
%   fval        :     MINUS value of the log posterior kernel at xparam1.
assia's avatar
assia committed
14 15 16 17 18 19 20 21
%   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
22
%
23

24
% Copyright (C) 2004-2009, 2010 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 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
%------------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
if options_.mode_compute ~= 1 & any(xparam1 < bayestopt_.lb)
    k = find(xparam1 < bayestopt_.lb);
    fval = bayestopt_.penalty+sum((bayestopt_.lb(k)-xparam1(k)).^2);
    cost_flag = 0;
    info = 41;
    return;
end
if options_.mode_compute ~= 1 & any(xparam1 > bayestopt_.ub)
    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
michel's avatar
michel committed
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
michel's avatar
michel committed
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);
michel's avatar
michel committed
75 76
    end
    offset = offset+estim_params_.nvn;
77 78
end
if estim_params_.ncx
michel's avatar
michel committed
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);
michel's avatar
michel committed
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
michel's avatar
michel committed
96 97
    end
    offset = offset+estim_params_.ncx;
98 99
end
if estim_params_.ncn 
michel's avatar
michel committed
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);
michel's avatar
michel committed
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
michel's avatar
michel committed
116 117
    end
    offset = offset+estim_params_.ncn;
118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149
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
%------------------------------------------------------------------------------
[T,R,SteadyState,info] = dynare_resolve(bayestopt_.restrict_var_list,...
                                        bayestopt_.restrict_columns,...
                                        bayestopt_.restrict_aux);
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)==6 ||info(1) == 19 || info(1) == 20 || info(1) == 21
    fval = bayestopt_.penalty+info(2);
    cost_flag = 0;
    return
end
bayestopt_.mf = bayestopt_.mf1;
if options_.noconstant
    constant = zeros(nobs,1);  
else    
    if options_.loglinear
        constant = log(SteadyState(bayestopt_.mfys));
    else
        constant = SteadyState(bayestopt_.mfys);
    end
end
if bayestopt_.with_trend
michel's avatar
michel committed
150
    trend_coeff = zeros(nobs,1);
151
    t = options_.trend_coeffs;
152
    for i=1:length(t)
153 154 155
        if ~isempty(t{i})
            trend_coeff(i) = evalin('base',t{i});
        end
michel's avatar
michel committed
156
    end
157
    trend = repmat(constant,1,gend)+trend_coeff*[1:gend];
158
else
159
    trend = repmat(constant,1,gend);
160 161 162 163 164 165 166 167 168
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;
169
if options_.lik_init == 1               % Kalman filter
170 171 172 173
    if kalman_algo ~= 2
        kalman_algo = 1;
    end
    Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);
174 175
    Pinf        = [];
elseif options_.lik_init == 2   % Old Diffuse Kalman filter
176 177 178 179 180
    if kalman_algo ~= 2
        kalman_algo = 1;
    end
    Pstar = options_.Harvey_scale_factor*eye(np);
    Pinf = [];
181
elseif options_.lik_init == 3   % Diffuse Kalman filter
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 309 310 311 312 313 314 315 316
    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; zeros(nk-size(Z,1),size(RR,2))])) < 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
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
    no_correlation_flag = 1;
    if length(H)==1 & H == 0
        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
    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
    no_correlation_flag = 1;
    if length(H)==1 & H == 0
        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
    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
317 318 319 320 321
if isnan(LIK)
    cost_flag = 0;
    return
end
if imag(LIK)~=0
322 323 324 325 326 327 328 329 330 331
    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);
options_.kalman_algo = kalman_algo;