From 035adeb89e5efe7d2dbe28bb312c53fdca39aacd Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Fri, 9 Oct 2015 14:20:11 +0200 Subject: [PATCH] Revert "More fixing related to objective_function_penalty_base" This reverts commit 1ad8df4635d99ba775e22f5693bf18c2cbb6c0aa. --- matlab/TaRB_optimizer_wrapper.m | 4 +- matlab/dsge_likelihood.m | 735 ++++++++++++++- matlab/dsge_likelihood_1.m | 858 ------------------ matlab/dsge_var_likelihood.m | 283 +++++- matlab/dsge_var_likelihood_1.m | 312 ------- matlab/dynare_estimation_1.m | 6 +- matlab/identification_analysis.m | 2 +- matlab/mode_check.m | 2 +- matlab/optimization/mr_hessian.m | 12 +- matlab/optimization/newrat.m | 6 +- .../optimization/penalty_objective_function.m | 6 +- matlab/osr1.m | 2 +- matlab/osr_obj.m | 60 +- matlab/osr_obj_1.m | 90 -- tests/estimation/TaRB/fs2000_tarb.mod | 1 - 15 files changed, 1084 insertions(+), 1295 deletions(-) delete mode 100644 matlab/dsge_likelihood_1.m delete mode 100644 matlab/dsge_var_likelihood_1.m delete mode 100644 matlab/osr_obj_1.m diff --git a/matlab/TaRB_optimizer_wrapper.m b/matlab/TaRB_optimizer_wrapper.m index ec91451b5..9c7ed4ccc 100644 --- a/matlab/TaRB_optimizer_wrapper.m +++ b/matlab/TaRB_optimizer_wrapper.m @@ -1,4 +1,4 @@ -function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff] = TaRB_optimizer_wrapper(optpar,par_vector,parameterindices,TargetFun,varargin) +function [fval,DLIK,Hess,exit_flag] = TaRB_optimizer_wrapper(optpar,par_vector,parameterindices,TargetFun,varargin) % function [fval,DLIK,Hess,exit_flag] = TaRB_optimizer_wrapper(optpar,par_vector,parameterindices,TargetFun,varargin) % Wrapper function for target function used in TaRB algorithm; reassembles % full parameter vector before calling target function @@ -36,5 +36,5 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff] = TaRB_optimize % along with Dynare. If not, see . par_vector(parameterindices,:)=optpar; %reassemble parameter -[fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff] = feval(TargetFun,par_vector,varargin{:}); %call target function +[fval,DLIK,Hess,exit_flag] = feval(TargetFun,par_vector,varargin{:}); %call target function diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m index 0552ae676..826c60496 100644 --- a/matlab/dsge_likelihood.m +++ b/matlab/dsge_likelihood.m @@ -1,5 +1,7 @@ -function [fval,DLIK,Hess,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,DynareResults] = dsge_likelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults,derivatives_info) -% Evaluates the posterior kernel of a dsge model. Deprecated interface. +function [fval,DLIK,Hess,exit_flag,SteadyState,trend_coeff,info,Model,DynareOptions,BayesInfo,DynareResults] = dsge_likelihood(xparam1,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults,derivatives_info) +% Evaluates the posterior kernel of a dsge model using the specified +% kalman_algo; the resulting posterior includes the 2*pi constant of the +% likelihood function %@info: %! @deftypefn {Function File} {[@var{fval},@var{exit_flag},@var{ys},@var{trend_coeff},@var{info},@var{Model},@var{DynareOptions},@var{BayesInfo},@var{DynareResults},@var{DLIK},@var{AHess}] =} dsge_likelihood (@var{xparam1},@var{DynareDataset},@var{DynareOptions},@var{Model},@var{EstimatedParameters},@var{BayesInfo},@var{DynareResults},@var{derivatives_flag}) @@ -37,7 +39,7 @@ function [fval,DLIK,Hess,exit_flag,ys,trend_coeff,info,Model,DynareOptions,Bayes %! Integer scalar, equal to zero if the routine return with a penalty (one otherwise). %! @item ys %! Vector of doubles, steady state level for the endogenous variables. -%! @item trend_coeffs +%! @item trend_coeff %! Matrix of doubles, coefficients of the deterministic trend in the measurement equation. %! @item info %! Integer scalar, error code. @@ -68,6 +70,8 @@ function [fval,DLIK,Hess,exit_flag,ys,trend_coeff,info,Model,DynareOptions,Bayes %! M_.params has been updated in the steadystate routine and has complex valued scalars. %! @item info==24 %! M_.params has been updated in the steadystate routine and has some NaNs. +%! @item info==26 +%! M_.params has been updated in the steadystate routine and has negative/0 values in loglinear model. %! @item info==30 %! Ergodic variance can't be computed. %! @item info==41 @@ -130,6 +134,725 @@ function [fval,DLIK,Hess,exit_flag,ys,trend_coeff,info,Model,DynareOptions,Bayes % AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT FR -[fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOptions,BayesInfo,DynareResults] = ... -dsge_likelihood_1(xparam1,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,... - BoundsInfo,DynareResults,derivatives_info); \ No newline at end of file +% Initialization of the returned variables and others... +fval = []; +SteadyState = []; +trend_coeff = []; +exit_flag = 1; +info = 0; +DLIK = []; +Hess = []; + +if DynareOptions.estimation_dll + [fval,exit_flag,SteadyState,trend_coeff,info,params,H,Q] ... + = logposterior(xparam1,DynareDataset, DynareOptions,Model, ... + EstimatedParameters,BayesInfo,DynareResults); + mexErrCheck('logposterior', exit_flag); + Model.params = params; + if ~isequal(Model.H,0) + Model.H = H; + end + Model.Sigma_e = Q; + DynareResults.dr.ys = SteadyState; + return +end + +% Set flag related to analytical derivatives. +analytic_derivation = DynareOptions.analytic_derivation; + +if analytic_derivation && DynareOptions.loglinear + error('The analytic_derivation and loglinear options are not compatible') +end + +if nargout==1, + analytic_derivation=0; +end + +if analytic_derivation, + kron_flag=DynareOptions.analytic_derivation_mode; +end + +%------------------------------------------------------------------------------ +% 1. Get the structural parameters & define penalties +%------------------------------------------------------------------------------ + +% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain. +if ~isequal(DynareOptions.mode_compute,1) && any(xparam1BoundsInfo.ub) + k = find(xparam1>BoundsInfo.ub); + fval = Inf; + exit_flag = 0; + info(1) = 42; + info(2) = sum((xparam1(k)-BoundsInfo.ub(k)).^2); + if analytic_derivation, + DLIK=ones(length(xparam1),1); + end + return +end + +% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H). +Model = set_all_parameters(xparam1,EstimatedParameters,Model); + +Q = Model.Sigma_e; +H = Model.H; + +% Test if Q is positive definite. +if ~issquare(Q) || EstimatedParameters.ncx || isfield(EstimatedParameters,'calibrated_covariances') + [Q_is_positive_definite, penalty] = ispd(Q); + if ~Q_is_positive_definite + fval = Inf; + exit_flag = 0; + info(1) = 43; + info(2) = penalty; + return + end + if isfield(EstimatedParameters,'calibrated_covariances') + correct_flag=check_consistency_covariances(Q); + if ~correct_flag + fval = Inf; + exit_flag = 0; + info(1) = 71; + info(2) = sum(Q(EstimatedParameters.calibrated_covariances.position).^2); + return + end + end +end + +% Test if H is positive definite. +if ~issquare(H) || EstimatedParameters.ncn || isfield(EstimatedParameters,'calibrated_covariances_ME') + [H_is_positive_definite, penalty] = ispd(H); + if ~H_is_positive_definite + fval = Inf; + exit_flag = 0; + info(1) = 44; + info(2) = penalty; + return + end + if isfield(EstimatedParameters,'calibrated_covariances_ME') + correct_flag=check_consistency_covariances(H); + if ~correct_flag + fval = Inf; + exit_flag = 0; + info(1) = 72; + info(2) = sum(H(EstimatedParameters.calibrated_covariances_ME.position).^2); + return + end + end +end + + +%------------------------------------------------------------------------------ +% 2. call model setup & reduction program +%------------------------------------------------------------------------------ + +% Linearize the model around the deterministic sdteadystate and extract the matrices of the state equation (T and R). +[T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict'); + +% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol). +if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) == 8 || ... + info(1) == 22 || info(1) == 24 || info(1) == 19 || info(1) == 25 || info(1) == 10 + fval = Inf; + info(2) = 0.1; + exit_flag = 0; + if analytic_derivation, + DLIK=ones(length(xparam1),1); + end + return +elseif info(1) == 3 || info(1) == 4 || info(1)==6 || info(1) == 20 || info(1) == 21 || info(1) == 23 || info(1)==26 + fval = Inf; + exit_flag = 0; + if analytic_derivation, + DLIK=ones(length(xparam1),1); + end + return +end + +% check endogenous prior restrictions +info=endogenous_prior_restrictions(T,R,Model,DynareOptions,DynareResults); +if info(1), + fval = Inf; + exit_flag = 0; + if analytic_derivation, + DLIK=ones(length(xparam1),1); + end + return +end +% + +% Define a vector of indices for the observed variables. Is this really usefull?... +BayesInfo.mf = BayesInfo.mf1; + +% Define the constant vector of the measurement equation. +if DynareOptions.noconstant + constant = zeros(DynareDataset.vobs,1); +else + if DynareOptions.loglinear + constant = log(SteadyState(BayesInfo.mfys)); + else + constant = SteadyState(BayesInfo.mfys); + end +end + +% Define the deterministic linear trend of the measurement equation. +if BayesInfo.with_trend + trend_coeff = zeros(DynareDataset.vobs,1); + t = DynareOptions.trend_coeffs; + for i=1:length(t) + if ~isempty(t{i}) + trend_coeff(i) = evalin('base',t{i}); + end + end + trend = repmat(constant,1,DynareDataset.nobs)+trend_coeff*[1:DynareDataset.nobs]; +else + trend_coeff = zeros(DynareDataset.vobs,1); + trend = repmat(constant,1,DynareDataset.nobs); +end + +% Get needed informations for kalman filter routines. +start = DynareOptions.presample+1; +Z = BayesInfo.mf; +no_missing_data_flag = ~DatasetInfo.missing.state; +mm = length(T); +pp = DynareDataset.vobs; +rr = length(Q); +kalman_tol = DynareOptions.kalman_tol; +diffuse_kalman_tol = DynareOptions.diffuse_kalman_tol; +riccati_tol = DynareOptions.riccati_tol; +Y = transpose(DynareDataset.data)-trend; + +%------------------------------------------------------------------------------ +% 3. Initial condition of the Kalman filter +%------------------------------------------------------------------------------ +kalman_algo = DynareOptions.kalman_algo; + +% resetting measurement errors covariance matrix for univariate filters +if (kalman_algo == 2) || (kalman_algo == 4) + if isequal(H,0) + H = zeros(pp,1); + mmm = mm; + else + if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal... + H = diag(H); + mmm = mm; + else + Z = [Z, eye(pp)]; + T = blkdiag(T,zeros(pp)); + Q = blkdiag(Q,H); + R = blkdiag(R,eye(pp)); + Pstar = blkdiag(Pstar,H); + Pinf = blkdiag(Pinf,zeros(pp)); + H = zeros(pp,1); + mmm = mm+pp; + end + end +end + + +diffuse_periods = 0; +correlated_errors_have_been_checked = 0; +singular_diffuse_filter = 0; +switch DynareOptions.lik_init + case 1% Standard initialization with the steady state of the state equation. + if kalman_algo~=2 + % Use standard kalman filter except if the univariate filter is explicitely choosen. + kalman_algo = 1; + end + if DynareOptions.lyapunov_fp == 1 + Pstar = lyapunov_symm(T,R*Q'*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 3, [], DynareOptions.debug); + elseif DynareOptions.lyapunov_db == 1 + Pstar = disclyap_fast(T,R*Q*R',DynareOptions.lyapunov_doubling_tol); + elseif DynareOptions.lyapunov_srs == 1 + Pstar = lyapunov_symm(T,Q,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 4, R, DynareOptions.debug); + else + Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); + end; + Pinf = []; + a = zeros(mm,1); + Zflag = 0; + case 2% Initialization with large numbers on the diagonal of the covariance matrix if the states (for non stationary models). + if kalman_algo ~= 2 + % Use standard kalman filter except if the univariate filter is explicitely choosen. + kalman_algo = 1; + end + Pstar = DynareOptions.Harvey_scale_factor*eye(mm); + Pinf = []; + a = zeros(mm,1); + Zflag = 0; + case 3% Diffuse Kalman filter (Durbin and Koopman) + % Use standard kalman filter except if the univariate filter is explicitely choosen. + if kalman_algo == 0 + kalman_algo = 3; + elseif ~((kalman_algo == 3) || (kalman_algo == 4)) + error(['The model requires Diffuse filter, but you specified a different Kalman filter. You must set options_.kalman_algo ' ... + 'to 0 (default), 3 or 4']) + end + [Z,T,R,QT,Pstar,Pinf] = schur_statespace_transformation(Z,T,R,Q,DynareOptions.qz_criterium); + Zflag = 1; + % Run diffuse kalman filter on first periods. + if (kalman_algo==3) + % Multivariate Diffuse Kalman Filter + if no_missing_data_flag + [dLIK,dlik,a,Pstar] = kalman_filter_d(Y, 1, size(Y,2), ... + zeros(mm,1), Pinf, Pstar, ... + kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ... + T,R,Q,H,Z,mm,pp,rr); + else + [dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations, ... + Y, 1, size(Y,2), ... + zeros(mm,1), Pinf, Pstar, ... + kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ... + T,R,Q,H,Z,mm,pp,rr); + end + diffuse_periods = length(dlik); + if isinf(dLIK) + % Go to univariate diffuse filter if singularity problem. + singular_diffuse_filter = 1; + end + end + if singular_diffuse_filter || (kalman_algo==4) + % Univariate Diffuse Kalman Filter + if isequal(H,0) + H1 = zeros(pp,1); + mmm = mm; + else + if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal... + H1 = diag(H); + mmm = mm; + else + Z = [Z, eye(pp)]; + T = blkdiag(T,zeros(pp)); + Q = blkdiag(Q,H); + R = blkdiag(R,eye(pp)); + Pstar = blkdiag(Pstar,H); + Pinf = blkdiag(Pinf,zeros(pp)); + H1 = zeros(pp,1); + mmm = mm+pp; + end + end + % no need to test again for correlation elements + correlated_errors_have_been_checked = 1; + + [dLIK,dlik,a,Pstar] = univariate_kalman_filter_d(DatasetInfo.missing.aindex,... + DatasetInfo.missing.number_of_observations,... + DatasetInfo.missing.no_more_missing_observations, ... + Y, 1, size(Y,2), ... + zeros(mmm,1), Pinf, Pstar, ... + kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ... + T,R,Q,H1,Z,mmm,pp,rr); + diffuse_periods = size(dlik,1); + end + if isnan(dLIK), + fval = Inf; + info(1) = 45; + info(2) = 0.1; + exit_flag = 0; + return + end + + case 4% Start from the solution of the Riccati equation. + if kalman_algo ~= 2 + kalman_algo = 1; + end + if isequal(H,0) + [err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,mm,length(Z)))); + else + [err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,mm,length(Z))),H); + end + if err + disp(['dsge_likelihood:: I am not able to solve the Riccati equation, so I switch to lik_init=1!']); + DynareOptions.lik_init = 1; + Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); + end + Pinf = []; + a = zeros(mm,1); + Zflag = 0; + case 5 % Old diffuse Kalman filter only for the non stationary variables + [eigenvect, eigenv] = eig(T); + eigenv = diag(eigenv); + nstable = length(find(abs(abs(eigenv)-1) > 1e-7)); + unstable = find(abs(abs(eigenv)-1) < 1e-7); + V = eigenvect(:,unstable); + indx_unstable = find(sum(abs(V),2)>1e-5); + stable = find(sum(abs(V),2)<1e-5); + nunit = length(eigenv) - nstable; + Pstar = options_.Harvey_scale_factor*eye(np); + if kalman_algo ~= 2 + kalman_algo = 1; + end + R_tmp = R(stable, :); + T_tmp = T(stable,stable); + if DynareOptions.lyapunov_fp == 1 + Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 3, [], DynareOptions.debug); + elseif DynareOptions.lyapunov_db == 1 + Pstar_tmp = disclyap_fast(T_tmp,R_tmp*Q*R_tmp',DynareOptions.lyapunov_doubling_tol); + elseif DynareOptions.lyapunov_srs == 1 + Pstar_tmp = lyapunov_symm(T_tmp,Q,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 4, R_tmp, DynareOptions.debug); + else + Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',DynareOptions.qz_criterium,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); + end + Pstar(stable, stable) = Pstar_tmp; + Pinf = []; + otherwise + error('dsge_likelihood:: Unknown initialization approach for the Kalman filter!') +end + +if analytic_derivation, + offset = EstimatedParameters.nvx; + offset = offset+EstimatedParameters.nvn; + offset = offset+EstimatedParameters.ncx; + offset = offset+EstimatedParameters.ncn; + + no_DLIK = 0; + full_Hess = analytic_derivation==2; + asy_Hess = analytic_derivation==-2; + outer_product_gradient = analytic_derivation==-1; + if asy_Hess, + analytic_derivation=1; + end + if outer_product_gradient, + analytic_derivation=1; + end + DLIK = []; + AHess = []; + iv = DynareResults.dr.restrict_var_list; + if nargin<10 || isempty(derivatives_info) + [A,B,nou,nou,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults); + if ~isempty(EstimatedParameters.var_exo) + indexo=EstimatedParameters.var_exo(:,1); + else + indexo=[]; + end + if ~isempty(EstimatedParameters.param_vals) + indparam=EstimatedParameters.param_vals(:,1); + else + indparam=[]; + end + + if full_Hess, + [dum, DT, DOm, DYss, dum2, D2T, D2Om, D2Yss] = getH(A, B, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,iv); + clear dum dum2; + else + [dum, DT, DOm, DYss] = getH(A, B, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,iv); + end + else + DT = derivatives_info.DT(iv,iv,:); + DOm = derivatives_info.DOm(iv,iv,:); + DYss = derivatives_info.DYss(iv,:); + if isfield(derivatives_info,'full_Hess'), + full_Hess = derivatives_info.full_Hess; + end + if full_Hess, + D2T = derivatives_info.D2T; + D2Om = derivatives_info.D2Om; + D2Yss = derivatives_info.D2Yss; + end + if isfield(derivatives_info,'no_DLIK'), + no_DLIK = derivatives_info.no_DLIK; + end + clear('derivatives_info'); + end + DYss = [zeros(size(DYss,1),offset) DYss]; + DH=zeros([length(H),length(H),length(xparam1)]); + DQ=zeros([size(Q),length(xparam1)]); + DP=zeros([size(T),length(xparam1)]); + if full_Hess, + for j=1:size(D2Yss,1), + tmp(j,:,:) = blkdiag(zeros(offset,offset), squeeze(D2Yss(j,:,:))); + end + D2Yss = tmp; + D2H=sparse(size(D2Om,1),size(D2Om,2)); %zeros([size(H),length(xparam1),length(xparam1)]); + D2P=sparse(size(D2Om,1),size(D2Om,2)); %zeros([size(T),length(xparam1),length(xparam1)]); + jcount=0; + end + if DynareOptions.lik_init==1, + for i=1:EstimatedParameters.nvx + k =EstimatedParameters.var_exo(i,1); + DQ(k,k,i) = 2*sqrt(Q(k,k)); + dum = lyapunov_symm(T,DOm(:,:,i),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug); +% kk = find(abs(dum) < 1e-12); +% dum(kk) = 0; + DP(:,:,i)=dum; + if full_Hess + for j=1:i, + jcount=jcount+1; + dum = lyapunov_symm(T,dyn_unvech(D2Om(:,jcount)),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug); +% kk = (abs(dum) < 1e-12); +% dum(kk) = 0; + D2P(:,jcount)=dyn_vech(dum); +% D2P(:,:,j,i)=dum; + end + end + end + end + offset = EstimatedParameters.nvx; + for i=1:EstimatedParameters.nvn + k = EstimatedParameters.var_endo(i,1); + DH(k,k,i+offset) = 2*sqrt(H(k,k)); + if full_Hess + D2H(k,k,i+offset,i+offset) = 2; + end + end + offset = offset + EstimatedParameters.nvn; + if DynareOptions.lik_init==1, + for j=1:EstimatedParameters.np + dum = lyapunov_symm(T,DT(:,:,j+offset)*Pstar*T'+T*Pstar*DT(:,:,j+offset)'+DOm(:,:,j+offset),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug); +% kk = find(abs(dum) < 1e-12); +% dum(kk) = 0; + DP(:,:,j+offset)=dum; + if full_Hess + DTj = DT(:,:,j+offset); + DPj = dum; + for i=1:j+offset, + jcount=jcount+1; + DTi = DT(:,:,i); + DPi = DP(:,:,i); + D2Tij = reshape(D2T(:,jcount),size(T)); + D2Omij = dyn_unvech(D2Om(:,jcount)); + tmp = D2Tij*Pstar*T' + T*Pstar*D2Tij' + DTi*DPj*T' + DTj*DPi*T' + T*DPj*DTi' + T*DPi*DTj' + DTi*Pstar*DTj' + DTj*Pstar*DTi' + D2Omij; + dum = lyapunov_symm(T,tmp,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug); +% dum(abs(dum)<1.e-12) = 0; + D2P(:,jcount) = dyn_vech(dum); +% D2P(:,:,j+offset,i) = dum; + end + end + end + end + if analytic_derivation==1, + analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP,asy_Hess}; + else + analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P}; + clear DT DYss DOm DP D2T D2Yss D2Om D2H D2P, + end +else + analytic_deriv_info={0}; +end + +%------------------------------------------------------------------------------ +% 4. Likelihood evaluation +%------------------------------------------------------------------------------ + +if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter + if no_missing_data_flag + if DynareOptions.block + [err, LIK] = block_kalman_filter(T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar); + mexErrCheck('block_kalman_filter', err); + else + [LIK,lik] = kalman_filter(Y,diffuse_periods+1,size(Y,2), ... + a,Pstar, ... + kalman_tol, riccati_tol, ... + DynareOptions.presample, ... + T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods, ... + analytic_deriv_info{:}); + + end + else + if 0 %DynareOptions.block + [err, LIK,lik] = block_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,... + T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar); + else + [LIK,lik] = missing_observations_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ... + a, Pstar, ... + kalman_tol, DynareOptions.riccati_tol, ... + DynareOptions.presample, ... + T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods); + end + end + if analytic_derivation, + LIK1=LIK; + LIK=LIK1{1}; + lik1=lik; + lik=lik1{1}; + end + if isinf(LIK) + if DynareOptions.use_univariate_filters_if_singularity_is_detected + if kalman_algo == 1 + kalman_algo = 2; + else + kalman_algo = 4; + end + else + if isinf(LIK) + fval = Inf; + info(1) = 66; + info(2) = 0.1; + exit_flag = 0; + return + end + end + else + if DynareOptions.lik_init==3 + LIK = LIK + dLIK; + if analytic_derivation==0 && nargout==2, + lik = [dlik; lik]; + end + end + end +end + +if (kalman_algo==2) || (kalman_algo==4) + % Univariate Kalman Filter + % resetting measurement error covariance matrix when necessary % + if ~correlated_errors_have_been_checked + if isequal(H,0) + H1 = zeros(pp,1); + mmm = mm; + if analytic_derivation, + DH = zeros(pp,length(xparam1)); + end + else + if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal... + H1 = diag(H); + mmm = mm; + clear tmp + if analytic_derivation, + for j=1:pp, + tmp(j,:)=DH(j,j,:); + end + DH=tmp; + end + else + Z = [Z, eye(pp)]; + T = blkdiag(T,zeros(pp)); + Q = blkdiag(Q,H); + R = blkdiag(R,eye(pp)); + Pstar = blkdiag(Pstar,H); + Pinf = blkdiag(Pinf,zeros(pp)); + H1 = zeros(pp,1); + mmm = mm+pp; + end + end + if analytic_derivation, + analytic_deriv_info{5}=DH; + end + end + + [LIK, lik] = univariate_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ... + a,Pstar, ... + DynareOptions.kalman_tol, ... + DynareOptions.riccati_tol, ... + DynareOptions.presample, ... + T,Q,R,H1,Z,mmm,pp,rr,Zflag,diffuse_periods,analytic_deriv_info{:}); + if analytic_derivation, + LIK1=LIK; + LIK=LIK1{1}; + lik1=lik; + lik=lik1{1}; + end + if DynareOptions.lik_init==3 + LIK = LIK+dLIK; + if analytic_derivation==0 && nargout==2, + lik = [dlik; lik]; + end + end +end + +if analytic_derivation + if no_DLIK==0 + DLIK = LIK1{2}; + % [DLIK] = score(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol); + end + if full_Hess , + Hess = -LIK1{3}; + % [Hess, DLL] = get_Hessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,Z,kalman_tol,riccati_tol); + % Hess0 = getHessian(Y,T,DT,D2T, R*Q*transpose(R),DOm,D2Om,Z,DYss,D2Yss); + end + if asy_Hess, +% if ~((kalman_algo==2) || (kalman_algo==4)), +% [Hess] = AHessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol); +% else + Hess = LIK1{3}; +% end + end +end + +if isnan(LIK) + fval = Inf; + info(1) = 45; + info(2) = 0.1; + exit_flag = 0; + return +end + +if imag(LIK)~=0 + fval = Inf; + info(1) = 46; + info(2) = 0.1; + exit_flag = 0; + return +end + +likelihood = LIK; + +% ------------------------------------------------------------------------------ +% 5. Adds prior if necessary +% ------------------------------------------------------------------------------ +if analytic_derivation + if full_Hess, + [lnprior, dlnprior, d2lnprior] = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4); + Hess = Hess - d2lnprior; + else + [lnprior, dlnprior] = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4); + end + if no_DLIK==0 + DLIK = DLIK - dlnprior'; + end + if outer_product_gradient, + dlik = lik1{2}; + dlik=[- dlnprior; dlik(start:end,:)]; + Hess = dlik'*dlik; + end +else + lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4); +end + +if DynareOptions.endogenous_prior==1 + if DynareOptions.lik_init==2 || DynareOptions.lik_init==3 + error('Endogenous prior not supported with non-stationary models') + else + [lnpriormom] = endogenous_prior(Y,Pstar,BayesInfo,H); + fval = (likelihood-lnprior-lnpriormom); + end +else + fval = (likelihood-lnprior); +end + +if DynareOptions.prior_restrictions.status + tmp = feval(DynareOptions.prior_restrictions.routine, Model, DynareResults, DynareOptions, DynareDataset, DatasetInfo); + fval = fval - tmp; +end + +if isnan(fval) + fval = Inf; + info(1) = 47; + info(2) = 0.1; + exit_flag = 0; + return +end + +if imag(fval)~=0 + fval = Inf; + info(1) = 48; + info(2) = 0.1; + exit_flag = 0; + return +end + +% Update DynareOptions.kalman_algo. +DynareOptions.kalman_algo = kalman_algo; + +if analytic_derivation==0 && nargout > 1 + lik=lik(start:end,:); + DLIK=[-lnprior; lik(:)]; +end \ No newline at end of file diff --git a/matlab/dsge_likelihood_1.m b/matlab/dsge_likelihood_1.m deleted file mode 100644 index cfa54c0e8..000000000 --- a/matlab/dsge_likelihood_1.m +++ /dev/null @@ -1,858 +0,0 @@ -function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOptions,BayesInfo,DynareResults] = dsge_likelihood_1(xparam1,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults,derivatives_info) -% Evaluates the posterior kernel of a dsge model using the specified -% kalman_algo; the resulting posterior includes the 2*pi constant of the -% likelihood function. - -%@info: -%! @deftypefn {Function File} {[@var{fval},@var{exit_flag},@var{ys},@var{trend_coeff},@var{info},@var{Model},@var{DynareOptions},@var{BayesInfo},@var{DynareResults},@var{DLIK},@var{AHess}] =} dsge_likelihood (@var{xparam1},@var{DynareDataset},@var{DynareOptions},@var{Model},@var{EstimatedParameters},@var{BayesInfo},@var{DynareResults},@var{derivatives_flag}) -%! @anchor{dsge_likelihood} -%! @sp 1 -%! Evaluates the posterior kernel of a dsge model. -%! @sp 2 -%! @strong{Inputs} -%! @sp 1 -%! @table @ @var -%! @item xparam1 -%! Vector of doubles, current values for the estimated parameters. -%! @item DynareDataset -%! Matlab's structure describing the dataset (initialized by dynare, see @ref{dataset_}). -%! @item DynareOptions -%! Matlab's structure describing the options (initialized by dynare, see @ref{options_}). -%! @item Model -%! Matlab's structure describing the Model (initialized by dynare, see @ref{M_}). -%! @item EstimatedParamemeters -%! Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}). -%! @item BayesInfo -%! Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}). -%! @item DynareResults -%! Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}). -%! @item derivates_flag -%! Integer scalar, flag for analytical derivatives of the likelihood. -%! @end table -%! @sp 2 -%! @strong{Outputs} -%! @sp 1 -%! @table @ @var -%! @item fval -%! Double scalar, value of (minus) the likelihood. -%! @item info -%! Integer scalar, error code. -%! @table @ @code -%! @item info==0 -%! No error. -%! @item info==1 -%! The model doesn't determine the current variables uniquely. -%! @item info==2 -%! MJDGGES returned an error code. -%! @item info==3 -%! Blanchard & Kahn conditions are not satisfied: no stable equilibrium. -%! @item info==4 -%! Blanchard & Kahn conditions are not satisfied: indeterminacy. -%! @item info==5 -%! Blanchard & Kahn conditions are not satisfied: indeterminacy due to rank failure. -%! @item info==6 -%! The jacobian evaluated at the deterministic steady state is complex. -%! @item info==19 -%! The steadystate routine thrown an exception (inconsistent deep parameters). -%! @item info==20 -%! Cannot find the steady state, info(2) contains the sum of square residuals (of the static equations). -%! @item info==21 -%! The steady state is complex, info(2) contains the sum of square of imaginary parts of the steady state. -%! @item info==22 -%! The steady has NaNs. -%! @item info==23 -%! M_.params has been updated in the steadystate routine and has complex valued scalars. -%! @item info==24 -%! M_.params has been updated in the steadystate routine and has some NaNs. -%! @item info==26 -%! M_.params has been updated in the steadystate routine and has negative/0 values in loglinear model. -%! @item info==30 -%! Ergodic variance can't be computed. -%! @item info==41 -%! At least one parameter is violating a lower bound condition. -%! @item info==42 -%! At least one parameter is violating an upper bound condition. -%! @item info==43 -%! The covariance matrix of the structural innovations is not positive definite. -%! @item info==44 -%! The covariance matrix of the measurement errors is not positive definite. -%! @item info==45 -%! Likelihood is not a number (NaN). -%! @item info==46 -%! Likelihood is a complex valued number. -%! @item info==47 -%! Posterior kernel is not a number (logged prior density is NaN) -%! @item info==48 -%! Posterior kernel is a complex valued number (logged prior density is complex). -%! @end table -%! @item exit_flag -%! Integer scalar, equal to zero if the routine return with a penalty (one otherwise). -%! @item ys -%! Vector of doubles, steady state level for the endogenous variables. -%! @item trend_coeff -%! Matrix of doubles, coefficients of the deterministic trend in the measurement equation. -%! @item Model -%! Matlab's structure describing the model (initialized by dynare, see @ref{M_}). -%! @item DynareOptions -%! Matlab's structure describing the options (initialized by dynare, see @ref{options_}). -%! @item BayesInfo -%! Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}). -%! @item DynareResults -%! Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}). -%! @item DLIK -%! Vector of doubles, score of the likelihood. -%! @item AHess -%! Matrix of doubles, asymptotic hessian matrix. -%! @end table -%! @sp 2 -%! @strong{This function is called by:} -%! @sp 1 -%! @ref{dynare_estimation_1}, @ref{mode_check} -%! @sp 2 -%! @strong{This function calls:} -%! @sp 1 -%! @ref{dynare_resolve}, @ref{lyapunov_symm}, @ref{schur_statespace_transformation}, @ref{kalman_filter_d}, @ref{missing_observations_kalman_filter_d}, @ref{univariate_kalman_filter_d}, @ref{kalman_steady_state}, @ref{getH}, @ref{kalman_filter}, @ref{score}, @ref{AHessian}, @ref{missing_observations_kalman_filter}, @ref{univariate_kalman_filter}, @ref{priordens} -%! @end deftypefn -%@eod: - -% Copyright (C) 2004-2013 Dynare Team -% -% 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 . - -% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT FR - -% Initialization of the returned variables and others... -fval = []; -SteadyState = []; -trend_coeff = []; -exit_flag = 1; -info = 0; -DLIK = []; -Hess = []; - -if DynareOptions.estimation_dll - [fval,exit_flag,SteadyState,trend_coeff,info,params,H,Q] ... - = logposterior(xparam1,DynareDataset, DynareOptions,Model, ... - EstimatedParameters,BayesInfo,DynareResults); - mexErrCheck('logposterior', exit_flag); - Model.params = params; - if ~isequal(Model.H,0) - Model.H = H; - end - Model.Sigma_e = Q; - DynareResults.dr.ys = SteadyState; - return -end - -% Set flag related to analytical derivatives. -analytic_derivation = DynareOptions.analytic_derivation; - -if analytic_derivation && DynareOptions.loglinear - error('The analytic_derivation and loglinear options are not compatible') -end - -if nargout==1, - analytic_derivation=0; -end - -if analytic_derivation, - kron_flag=DynareOptions.analytic_derivation_mode; -end - -%------------------------------------------------------------------------------ -% 1. Get the structural parameters & define penalties -%------------------------------------------------------------------------------ - -% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain. -if ~isequal(DynareOptions.mode_compute,1) && any(xparam1BoundsInfo.ub) - k = find(xparam1>BoundsInfo.ub); - fval = Inf; - exit_flag = 0; - info(1) = 42; - info(2) = sum((xparam1(k)-BoundsInfo.ub(k)).^2); - if analytic_derivation, - DLIK=ones(length(xparam1),1); - end - return -end - -% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H). -Model = set_all_parameters(xparam1,EstimatedParameters,Model); - -Q = Model.Sigma_e; -H = Model.H; - -% Test if Q is positive definite. -if ~issquare(Q) || EstimatedParameters.ncx || isfield(EstimatedParameters,'calibrated_covariances') - [Q_is_positive_definite, penalty] = ispd(Q); - if ~Q_is_positive_definite - fval = Inf; - exit_flag = 0; - info(1) = 43; - info(2) = penalty; - return - end - if isfield(EstimatedParameters,'calibrated_covariances') - correct_flag=check_consistency_covariances(Q); - if ~correct_flag - fval = Inf; - exit_flag = 0; - info(1) = 71; - info(2) = sum(Q(EstimatedParameters.calibrated_covariances.position).^2); - return - end - end -end - -% Test if H is positive definite. -if ~issquare(H) || EstimatedParameters.ncn || isfield(EstimatedParameters,'calibrated_covariances_ME') - [H_is_positive_definite, penalty] = ispd(H); - if ~H_is_positive_definite - fval = Inf; - exit_flag = 0; - info(1) = 44; - info(2) = penalty; - return - end - if isfield(EstimatedParameters,'calibrated_covariances_ME') - correct_flag=check_consistency_covariances(H); - if ~correct_flag - fval = Inf; - exit_flag = 0; - info(1) = 72; - info(2) = sum(H(EstimatedParameters.calibrated_covariances_ME.position).^2); - return - end - end -end - - -%------------------------------------------------------------------------------ -% 2. call model setup & reduction program -%------------------------------------------------------------------------------ - -% Linearize the model around the deterministic sdteadystate and extract the matrices of the state equation (T and R). -[T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict'); - -% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol). -if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) == 8 || ... - info(1) == 22 || info(1) == 24 || info(1) == 19 || info(1) == 25 || info(1) == 10 - fval = Inf; - info(2) = 0.1; - exit_flag = 0; - if analytic_derivation, - DLIK=ones(length(xparam1),1); - end - return -elseif info(1) == 3 || info(1) == 4 || info(1)==6 || info(1) == 20 || info(1) == 21 || info(1) == 23 || info(1)==26 - fval = Inf; - exit_flag = 0; - if analytic_derivation, - DLIK=ones(length(xparam1),1); - end - return -end - -% check endogenous prior restrictions -info=endogenous_prior_restrictions(T,R,Model,DynareOptions,DynareResults); -if info(1), - fval = Inf; - exit_flag = 0; - if analytic_derivation, - DLIK=ones(length(xparam1),1); - end - return -end -% - -% Define a vector of indices for the observed variables. Is this really usefull?... -BayesInfo.mf = BayesInfo.mf1; - -% Define the constant vector of the measurement equation. -if DynareOptions.noconstant - constant = zeros(DynareDataset.vobs,1); -else - if DynareOptions.loglinear - constant = log(SteadyState(BayesInfo.mfys)); - else - constant = SteadyState(BayesInfo.mfys); - end -end - -% Define the deterministic linear trend of the measurement equation. -if BayesInfo.with_trend - trend_coeff = zeros(DynareDataset.vobs,1); - t = DynareOptions.trend_coeffs; - for i=1:length(t) - if ~isempty(t{i}) - trend_coeff(i) = evalin('base',t{i}); - end - end - trend = repmat(constant,1,DynareDataset.nobs)+trend_coeff*[1:DynareDataset.nobs]; -else - trend_coeff = zeros(DynareDataset.vobs,1); - trend = repmat(constant,1,DynareDataset.nobs); -end - -% Get needed informations for kalman filter routines. -start = DynareOptions.presample+1; -Z = BayesInfo.mf; -no_missing_data_flag = ~DatasetInfo.missing.state; -mm = length(T); -pp = DynareDataset.vobs; -rr = length(Q); -kalman_tol = DynareOptions.kalman_tol; -diffuse_kalman_tol = DynareOptions.diffuse_kalman_tol; -riccati_tol = DynareOptions.riccati_tol; -Y = transpose(DynareDataset.data)-trend; - -%------------------------------------------------------------------------------ -% 3. Initial condition of the Kalman filter -%------------------------------------------------------------------------------ -kalman_algo = DynareOptions.kalman_algo; - -% resetting measurement errors covariance matrix for univariate filters -if (kalman_algo == 2) || (kalman_algo == 4) - if isequal(H,0) - H = zeros(pp,1); - mmm = mm; - else - if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal... - H = diag(H); - mmm = mm; - else - Z = [Z, eye(pp)]; - T = blkdiag(T,zeros(pp)); - Q = blkdiag(Q,H); - R = blkdiag(R,eye(pp)); - Pstar = blkdiag(Pstar,H); - Pinf = blkdiag(Pinf,zeros(pp)); - H = zeros(pp,1); - mmm = mm+pp; - end - end -end - - -diffuse_periods = 0; -correlated_errors_have_been_checked = 0; -singular_diffuse_filter = 0; -switch DynareOptions.lik_init - case 1% Standard initialization with the steady state of the state equation. - if kalman_algo~=2 - % Use standard kalman filter except if the univariate filter is explicitely choosen. - kalman_algo = 1; - end - if DynareOptions.lyapunov_fp == 1 - Pstar = lyapunov_symm(T,R*Q'*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 3, [], DynareOptions.debug); - elseif DynareOptions.lyapunov_db == 1 - Pstar = disclyap_fast(T,R*Q*R',DynareOptions.lyapunov_doubling_tol); - elseif DynareOptions.lyapunov_srs == 1 - Pstar = lyapunov_symm(T,Q,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 4, R, DynareOptions.debug); - else - Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); - end; - Pinf = []; - a = zeros(mm,1); - Zflag = 0; - case 2% Initialization with large numbers on the diagonal of the covariance matrix if the states (for non stationary models). - if kalman_algo ~= 2 - % Use standard kalman filter except if the univariate filter is explicitely choosen. - kalman_algo = 1; - end - Pstar = DynareOptions.Harvey_scale_factor*eye(mm); - Pinf = []; - a = zeros(mm,1); - Zflag = 0; - case 3% Diffuse Kalman filter (Durbin and Koopman) - % Use standard kalman filter except if the univariate filter is explicitely choosen. - if kalman_algo == 0 - kalman_algo = 3; - elseif ~((kalman_algo == 3) || (kalman_algo == 4)) - error(['The model requires Diffuse filter, but you specified a different Kalman filter. You must set options_.kalman_algo ' ... - 'to 0 (default), 3 or 4']) - end - [Z,T,R,QT,Pstar,Pinf] = schur_statespace_transformation(Z,T,R,Q,DynareOptions.qz_criterium); - Zflag = 1; - % Run diffuse kalman filter on first periods. - if (kalman_algo==3) - % Multivariate Diffuse Kalman Filter - if no_missing_data_flag - [dLIK,dlik,a,Pstar] = kalman_filter_d(Y, 1, size(Y,2), ... - zeros(mm,1), Pinf, Pstar, ... - kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ... - T,R,Q,H,Z,mm,pp,rr); - else - [dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations, ... - Y, 1, size(Y,2), ... - zeros(mm,1), Pinf, Pstar, ... - kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ... - T,R,Q,H,Z,mm,pp,rr); - end - diffuse_periods = length(dlik); - if isinf(dLIK) - % Go to univariate diffuse filter if singularity problem. - singular_diffuse_filter = 1; - end - end - if singular_diffuse_filter || (kalman_algo==4) - % Univariate Diffuse Kalman Filter - if isequal(H,0) - H1 = zeros(pp,1); - mmm = mm; - else - if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal... - H1 = diag(H); - mmm = mm; - else - Z = [Z, eye(pp)]; - T = blkdiag(T,zeros(pp)); - Q = blkdiag(Q,H); - R = blkdiag(R,eye(pp)); - Pstar = blkdiag(Pstar,H); - Pinf = blkdiag(Pinf,zeros(pp)); - H1 = zeros(pp,1); - mmm = mm+pp; - end - end - % no need to test again for correlation elements - correlated_errors_have_been_checked = 1; - - [dLIK,dlik,a,Pstar] = univariate_kalman_filter_d(DatasetInfo.missing.aindex,... - DatasetInfo.missing.number_of_observations,... - DatasetInfo.missing.no_more_missing_observations, ... - Y, 1, size(Y,2), ... - zeros(mmm,1), Pinf, Pstar, ... - kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ... - T,R,Q,H1,Z,mmm,pp,rr); - diffuse_periods = size(dlik,1); - end - if isnan(dLIK), - fval = Inf; - info(1) = 45; - info(2) = 0.1; - exit_flag = 0; - return - end - - case 4% Start from the solution of the Riccati equation. - if kalman_algo ~= 2 - kalman_algo = 1; - end - if isequal(H,0) - [err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,mm,length(Z)))); - else - [err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,mm,length(Z))),H); - end - if err - disp(['dsge_likelihood:: I am not able to solve the Riccati equation, so I switch to lik_init=1!']); - DynareOptions.lik_init = 1; - Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); - end - Pinf = []; - a = zeros(mm,1); - Zflag = 0; - case 5 % Old diffuse Kalman filter only for the non stationary variables - [eigenvect, eigenv] = eig(T); - eigenv = diag(eigenv); - nstable = length(find(abs(abs(eigenv)-1) > 1e-7)); - unstable = find(abs(abs(eigenv)-1) < 1e-7); - V = eigenvect(:,unstable); - indx_unstable = find(sum(abs(V),2)>1e-5); - stable = find(sum(abs(V),2)<1e-5); - nunit = length(eigenv) - nstable; - Pstar = options_.Harvey_scale_factor*eye(np); - if kalman_algo ~= 2 - kalman_algo = 1; - end - R_tmp = R(stable, :); - T_tmp = T(stable,stable); - if DynareOptions.lyapunov_fp == 1 - Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 3, [], DynareOptions.debug); - elseif DynareOptions.lyapunov_db == 1 - Pstar_tmp = disclyap_fast(T_tmp,R_tmp*Q*R_tmp',DynareOptions.lyapunov_doubling_tol); - elseif DynareOptions.lyapunov_srs == 1 - Pstar_tmp = lyapunov_symm(T_tmp,Q,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 4, R_tmp, DynareOptions.debug); - else - Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',DynareOptions.qz_criterium,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); - end - Pstar(stable, stable) = Pstar_tmp; - Pinf = []; - otherwise - error('dsge_likelihood:: Unknown initialization approach for the Kalman filter!') -end - -if analytic_derivation, - offset = EstimatedParameters.nvx; - offset = offset+EstimatedParameters.nvn; - offset = offset+EstimatedParameters.ncx; - offset = offset+EstimatedParameters.ncn; - - no_DLIK = 0; - full_Hess = analytic_derivation==2; - asy_Hess = analytic_derivation==-2; - outer_product_gradient = analytic_derivation==-1; - if asy_Hess, - analytic_derivation=1; - end - if outer_product_gradient, - analytic_derivation=1; - end - DLIK = []; - AHess = []; - iv = DynareResults.dr.restrict_var_list; - if nargin<10 || isempty(derivatives_info) - [A,B,nou,nou,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults); - if ~isempty(EstimatedParameters.var_exo) - indexo=EstimatedParameters.var_exo(:,1); - else - indexo=[]; - end - if ~isempty(EstimatedParameters.param_vals) - indparam=EstimatedParameters.param_vals(:,1); - else - indparam=[]; - end - - if full_Hess, - [dum, DT, DOm, DYss, dum2, D2T, D2Om, D2Yss] = getH(A, B, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,iv); - clear dum dum2; - else - [dum, DT, DOm, DYss] = getH(A, B, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,iv); - end - else - DT = derivatives_info.DT(iv,iv,:); - DOm = derivatives_info.DOm(iv,iv,:); - DYss = derivatives_info.DYss(iv,:); - if isfield(derivatives_info,'full_Hess'), - full_Hess = derivatives_info.full_Hess; - end - if full_Hess, - D2T = derivatives_info.D2T; - D2Om = derivatives_info.D2Om; - D2Yss = derivatives_info.D2Yss; - end - if isfield(derivatives_info,'no_DLIK'), - no_DLIK = derivatives_info.no_DLIK; - end - clear('derivatives_info'); - end - DYss = [zeros(size(DYss,1),offset) DYss]; - DH=zeros([length(H),length(H),length(xparam1)]); - DQ=zeros([size(Q),length(xparam1)]); - DP=zeros([size(T),length(xparam1)]); - if full_Hess, - for j=1:size(D2Yss,1), - tmp(j,:,:) = blkdiag(zeros(offset,offset), squeeze(D2Yss(j,:,:))); - end - D2Yss = tmp; - D2H=sparse(size(D2Om,1),size(D2Om,2)); %zeros([size(H),length(xparam1),length(xparam1)]); - D2P=sparse(size(D2Om,1),size(D2Om,2)); %zeros([size(T),length(xparam1),length(xparam1)]); - jcount=0; - end - if DynareOptions.lik_init==1, - for i=1:EstimatedParameters.nvx - k =EstimatedParameters.var_exo(i,1); - DQ(k,k,i) = 2*sqrt(Q(k,k)); - dum = lyapunov_symm(T,DOm(:,:,i),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug); -% kk = find(abs(dum) < 1e-12); -% dum(kk) = 0; - DP(:,:,i)=dum; - if full_Hess - for j=1:i, - jcount=jcount+1; - dum = lyapunov_symm(T,dyn_unvech(D2Om(:,jcount)),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug); -% kk = (abs(dum) < 1e-12); -% dum(kk) = 0; - D2P(:,jcount)=dyn_vech(dum); -% D2P(:,:,j,i)=dum; - end - end - end - end - offset = EstimatedParameters.nvx; - for i=1:EstimatedParameters.nvn - k = EstimatedParameters.var_endo(i,1); - DH(k,k,i+offset) = 2*sqrt(H(k,k)); - if full_Hess - D2H(k,k,i+offset,i+offset) = 2; - end - end - offset = offset + EstimatedParameters.nvn; - if DynareOptions.lik_init==1, - for j=1:EstimatedParameters.np - dum = lyapunov_symm(T,DT(:,:,j+offset)*Pstar*T'+T*Pstar*DT(:,:,j+offset)'+DOm(:,:,j+offset),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug); -% kk = find(abs(dum) < 1e-12); -% dum(kk) = 0; - DP(:,:,j+offset)=dum; - if full_Hess - DTj = DT(:,:,j+offset); - DPj = dum; - for i=1:j+offset, - jcount=jcount+1; - DTi = DT(:,:,i); - DPi = DP(:,:,i); - D2Tij = reshape(D2T(:,jcount),size(T)); - D2Omij = dyn_unvech(D2Om(:,jcount)); - tmp = D2Tij*Pstar*T' + T*Pstar*D2Tij' + DTi*DPj*T' + DTj*DPi*T' + T*DPj*DTi' + T*DPi*DTj' + DTi*Pstar*DTj' + DTj*Pstar*DTi' + D2Omij; - dum = lyapunov_symm(T,tmp,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug); -% dum(abs(dum)<1.e-12) = 0; - D2P(:,jcount) = dyn_vech(dum); -% D2P(:,:,j+offset,i) = dum; - end - end - end - end - if analytic_derivation==1, - analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP,asy_Hess}; - else - analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P}; - clear DT DYss DOm DP D2T D2Yss D2Om D2H D2P, - end -else - analytic_deriv_info={0}; -end - -%------------------------------------------------------------------------------ -% 4. Likelihood evaluation -%------------------------------------------------------------------------------ - -if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter - if no_missing_data_flag - if DynareOptions.block - [err, LIK] = block_kalman_filter(T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar); - mexErrCheck('block_kalman_filter', err); - else - [LIK,lik] = kalman_filter(Y,diffuse_periods+1,size(Y,2), ... - a,Pstar, ... - kalman_tol, riccati_tol, ... - DynareOptions.presample, ... - T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods, ... - analytic_deriv_info{:}); - - end - else - if 0 %DynareOptions.block - [err, LIK,lik] = block_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,... - T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar); - else - [LIK,lik] = missing_observations_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ... - a, Pstar, ... - kalman_tol, DynareOptions.riccati_tol, ... - DynareOptions.presample, ... - T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods); - end - end - if analytic_derivation, - LIK1=LIK; - LIK=LIK1{1}; - lik1=lik; - lik=lik1{1}; - end - if isinf(LIK) - if DynareOptions.use_univariate_filters_if_singularity_is_detected - if kalman_algo == 1 - kalman_algo = 2; - else - kalman_algo = 4; - end - else - if isinf(LIK) - fval = Inf; - info(1) = 66; - info(2) = 0.1; - exit_flag = 0; - return - end - end - else - if DynareOptions.lik_init==3 - LIK = LIK + dLIK; - if analytic_derivation==0 && nargout==2, - lik = [dlik; lik]; - end - end - end -end - -if (kalman_algo==2) || (kalman_algo==4) - % Univariate Kalman Filter - % resetting measurement error covariance matrix when necessary % - if ~correlated_errors_have_been_checked - if isequal(H,0) - H1 = zeros(pp,1); - mmm = mm; - if analytic_derivation, - DH = zeros(pp,length(xparam1)); - end - else - if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal... - H1 = diag(H); - mmm = mm; - clear tmp - if analytic_derivation, - for j=1:pp, - tmp(j,:)=DH(j,j,:); - end - DH=tmp; - end - else - Z = [Z, eye(pp)]; - T = blkdiag(T,zeros(pp)); - Q = blkdiag(Q,H); - R = blkdiag(R,eye(pp)); - Pstar = blkdiag(Pstar,H); - Pinf = blkdiag(Pinf,zeros(pp)); - H1 = zeros(pp,1); - mmm = mm+pp; - end - end - if analytic_derivation, - analytic_deriv_info{5}=DH; - end - end - - [LIK, lik] = univariate_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ... - a,Pstar, ... - DynareOptions.kalman_tol, ... - DynareOptions.riccati_tol, ... - DynareOptions.presample, ... - T,Q,R,H1,Z,mmm,pp,rr,Zflag,diffuse_periods,analytic_deriv_info{:}); - if analytic_derivation, - LIK1=LIK; - LIK=LIK1{1}; - lik1=lik; - lik=lik1{1}; - end - if DynareOptions.lik_init==3 - LIK = LIK+dLIK; - if analytic_derivation==0 && nargout==2, - lik = [dlik; lik]; - end - end -end - -if analytic_derivation - if no_DLIK==0 - DLIK = LIK1{2}; - % [DLIK] = score(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol); - end - if full_Hess , - Hess = -LIK1{3}; - % [Hess, DLL] = get_Hessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,Z,kalman_tol,riccati_tol); - % Hess0 = getHessian(Y,T,DT,D2T, R*Q*transpose(R),DOm,D2Om,Z,DYss,D2Yss); - end - if asy_Hess, -% if ~((kalman_algo==2) || (kalman_algo==4)), -% [Hess] = AHessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol); -% else - Hess = LIK1{3}; -% end - end -end - -if isnan(LIK) - fval = Inf; - info(1) = 45; - info(2) = 0.1; - exit_flag = 0; - return -end - -if imag(LIK)~=0 - fval = Inf; - info(1) = 46; - info(2) = 0.1; - exit_flag = 0; - return -end - -likelihood = LIK; - -% ------------------------------------------------------------------------------ -% 5. Adds prior if necessary -% ------------------------------------------------------------------------------ -if analytic_derivation - if full_Hess, - [lnprior, dlnprior, d2lnprior] = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4); - Hess = Hess - d2lnprior; - else - [lnprior, dlnprior] = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4); - end - if no_DLIK==0 - DLIK = DLIK - dlnprior'; - end - if outer_product_gradient, - dlik = lik1{2}; - dlik=[- dlnprior; dlik(start:end,:)]; - Hess = dlik'*dlik; - end -else - lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4); -end - -if DynareOptions.endogenous_prior==1 - if DynareOptions.lik_init==2 || DynareOptions.lik_init==3 - error('Endogenous prior not supported with non-stationary models') - else - [lnpriormom] = endogenous_prior(Y,Pstar,BayesInfo,H); - fval = (likelihood-lnprior-lnpriormom); - end -else - fval = (likelihood-lnprior); -end - -if DynareOptions.prior_restrictions.status - tmp = feval(DynareOptions.prior_restrictions.routine, Model, DynareResults, DynareOptions, DynareDataset, DatasetInfo); - fval = fval - tmp; -end - -if isnan(fval) - fval = Inf; - info(1) = 47; - info(2) = 0.1; - exit_flag = 0; - return -end - -if imag(fval)~=0 - fval = Inf; - info(1) = 48; - info(2) = 0.1; - exit_flag = 0; - return -end - -% Update DynareOptions.kalman_algo. -DynareOptions.kalman_algo = kalman_algo; - -if analytic_derivation==0 && nargout > 1 - lik=lik(start:end,:); - DLIK=[-lnprior; lik(:)]; -end \ No newline at end of file diff --git a/matlab/dsge_var_likelihood.m b/matlab/dsge_var_likelihood.m index 044ed7269..d0a5a8b4d 100644 --- a/matlab/dsge_var_likelihood.m +++ b/matlab/dsge_var_likelihood.m @@ -1,5 +1,5 @@ -function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) -% Evaluates the posterior kernel of the bvar-dsge model. Deprecated interface. +function [fval,grad,hess,exit_flag,SteadyState,trend_coeff,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(xparam1,DynareDataset,DynareInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults) +% Evaluates the posterior kernel of the bvar-dsge model. % % INPUTS % o xparam1 [double] Vector of model's parameters. @@ -8,6 +8,10 @@ function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelih % 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 SteadyState [double] Steady state vector possibly recomputed +% by call to dynare_results() +% o trend_coeff [double] place holder for trend coefficients, +% currently not supported by dsge_var % 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). @@ -34,6 +38,275 @@ function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelih % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -[fval,info,exit_flag,grad,hess,SteadyState,trend_coeff,PHI,SIGMAu,iXX,prior] = ... - dsge_var_likelihood_1(xparam1,DynareDataset,DynareInfo,DynareOptions,Model,... - EstimatedParameters,BayesInfo,BoundsInfo,DynareResults); \ No newline at end of file +persistent dsge_prior_weight_idx + +grad=[]; +hess=[]; +exit_flag = []; +info = 0; +PHI = []; +SIGMAu = []; +iXX = []; +prior = []; +SteadyState = []; +trend_coeff = []; + +% Initialization of of the index for parameter dsge_prior_weight in Model.params. +if isempty(dsge_prior_weight_idx) + dsge_prior_weight_idx = strmatch('dsge_prior_weight',Model.param_names); +end + +% Get the number of estimated (dsge) parameters. +nx = EstimatedParameters.nvx + EstimatedParameters.np; + +% Get the number of observed variables in the VAR model. +NumberOfObservedVariables = DynareDataset.vobs; + +% Get the number of observations. +NumberOfObservations = DynareDataset.nobs; + + +% Get the number of lags in the VAR model. +NumberOfLags = DynareOptions.dsge_varlag; + +% Get the number of parameters in the VAR model. +NumberOfParameters = NumberOfObservedVariables*NumberOfLags ; +if ~DynareOptions.noconstant + NumberOfParameters = NumberOfParameters + 1; +end + +% Get empirical second order moments for the observed variables. +mYY = evalin('base', 'mYY'); +mYX = evalin('base', 'mYX'); +mXY = evalin('base', 'mXY'); +mXX = evalin('base', 'mXX'); + +% Initialize some of the output arguments. +fval = []; +exit_flag = 1; + +% Return, with endogenous penalty, if some dsge-parameters are smaller than the lower bound of the prior domain. +if DynareOptions.mode_compute ~= 1 && any(xparam1 < BoundsInfo.lb) + fval = Inf; + exit_flag = 0; + info(1) = 41; + k = find(xparam1 < BoundsInfo.lb); + info(2) = sum((BoundsInfo.lb(k)-xparam1(k)).^2); + return; +end + +% Return, with endogenous penalty, if some dsge-parameters are greater than the upper bound of the prior domain. +if DynareOptions.mode_compute ~= 1 && any(xparam1 > BoundsInfo.ub) + fval = Inf; + exit_flag = 0; + info(1) = 42; + k = find(xparam1 > BoundsInfo.ub); + info(2) = sum((xparam1(k)-BoundsInfo.ub(k)).^2); + return; +end + +% Get the variance of each structural innovation. +Q = Model.Sigma_e; +for i=1:EstimatedParameters.nvx + k = EstimatedParameters.var_exo(i,1); + Q(k,k) = xparam1(i)*xparam1(i); +end +offset = EstimatedParameters.nvx; + +% Update Model.params and Model.Sigma_e. +Model.params(EstimatedParameters.param_vals(:,1)) = xparam1(offset+1:end); +Model.Sigma_e = Q; + +% Get the weight of the dsge prior. +dsge_prior_weight = Model.params(dsge_prior_weight_idx); + +% Is the dsge prior proper? +if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/ ... + NumberOfObservations; + fval = Inf; + exit_flag = 0; + info(1) = 51; + info(2) = abs(NumberOfObservations*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables)); + % info(2)=dsge_prior_weight; + % info(3)=(NumberOfParameters+NumberOfObservedVariables)/DynareDataset.nobs; + return +end + +%------------------------------------------------------------------------------ +% 2. call model setup & reduction program +%------------------------------------------------------------------------------ + +% Solve the Dsge model and get the matrices of the reduced form solution. T and R are the matrices of the +% state equation +[T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict'); + +% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol). +if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) == 8 || ... + info(1) == 22 || info(1) == 24 || info(1) == 25 || info(1) == 10 + fval = Inf; + info(2) = 0.1; + exit_flag = 0; + return +elseif info(1) == 3 || info(1) == 4 || info(1) == 19 || info(1) == 20 || info(1) == 21 + fval = Inf; + info(2) = 0.1; + exit_flag = 0; + return +end + +% Define the mean/steady state vector. +if ~DynareOptions.noconstant + if DynareOptions.loglinear + constant = transpose(log(SteadyState(BayesInfo.mfys))); + else + constant = transpose(SteadyState(BayesInfo.mfys)); + end +else + constant = zeros(1,NumberOfObservedVariables); +end + + +%------------------------------------------------------------------------------ +% 3. theoretical moments (second order) +%------------------------------------------------------------------------------ + +% Compute the theoretical second order moments +tmp0 = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); +mf = BayesInfo.mf1; + +% Get the non centered second order moments +TheoreticalAutoCovarianceOfTheObservedVariables = zeros(NumberOfObservedVariables,NumberOfObservedVariables,NumberOfLags+1); +TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) = tmp0(mf,mf)+constant'*constant; +for lag = 1:NumberOfLags + tmp0 = T*tmp0; + TheoreticalAutoCovarianceOfTheObservedVariables(:,:,lag+1) = tmp0(mf,mf) + constant'*constant; +end + +% Build the theoretical "covariance" between Y and X +GYX = zeros(NumberOfObservedVariables,NumberOfParameters); +for i=1:NumberOfLags + GYX(:,(i-1)*NumberOfObservedVariables+1:i*NumberOfObservedVariables) = TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1); +end +if ~DynareOptions.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 ~DynareOptions.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)% Evaluation of the likelihood of the dsge-var model when the dsge prior weight is finite. + tmp0 = dsge_prior_weight*NumberOfObservations*TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) + mYY ; + tmp1 = dsge_prior_weight*NumberOfObservations*GYX + mYX; + tmp2 = inv(dsge_prior_weight*NumberOfObservations*GXX+mXX); + SIGMAu = tmp0 - tmp1*tmp2*tmp1'; clear('tmp0'); + [SIGMAu_is_positive_definite, penalty] = ispd(SIGMAu); + if ~SIGMAu_is_positive_definite + fval = Inf; + info(1) = 52; + info(2) = penalty; + exit_flag = 0; + return; + end + SIGMAu = SIGMAu / (NumberOfObservations*(1+dsge_prior_weight)); + PHI = tmp2*tmp1'; clear('tmp1'); + prodlng1 = sum(gammaln(.5*((1+dsge_prior_weight)*NumberOfObservations- ... + NumberOfObservedVariables*NumberOfLags ... + +1-(1:NumberOfObservedVariables)'))); + prodlng2 = sum(gammaln(.5*(dsge_prior_weight*NumberOfObservations- ... + NumberOfObservedVariables*NumberOfLags ... + +1-(1:NumberOfObservedVariables)'))); + lik = .5*NumberOfObservedVariables*log(det(dsge_prior_weight*NumberOfObservations*GXX+mXX)) ... + + .5*((dsge_prior_weight+1)*NumberOfObservations-NumberOfParameters)*log(det((dsge_prior_weight+1)*NumberOfObservations*SIGMAu)) ... + - .5*NumberOfObservedVariables*log(det(dsge_prior_weight*NumberOfObservations*GXX)) ... + - .5*(dsge_prior_weight*NumberOfObservations-NumberOfParameters)*log(det(dsge_prior_weight*NumberOfObservations*(GYY-GYX*inv(GXX)*GYX'))) ... + + .5*NumberOfObservedVariables*NumberOfObservations*log(2*pi) ... + - .5*log(2)*NumberOfObservedVariables*((dsge_prior_weight+1)*NumberOfObservations-NumberOfParameters) ... + + .5*log(2)*NumberOfObservedVariables*(dsge_prior_weight*NumberOfObservations-NumberOfParameters) ... + - prodlng1 + prodlng2; +else% Evaluation of the likelihood of the dsge-var model when the dsge prior weight is infinite. + iGXX = inv(GXX); + SIGMAu = GYY - GYX*iGXX*transpose(GYX); + PHI = iGXX*transpose(GYX); + lik = NumberOfObservations * ( log(det(SIGMAu)) + NumberOfObservedVariables*log(2*pi) + ... + trace(inv(SIGMAu)*(mYY - transpose(mYX*PHI) - mYX*PHI + transpose(PHI)*mXX*PHI)/NumberOfObservations)); + lik = .5*lik;% Minus likelihood +end + +if isnan(lik) + info(1) = 45; + info(2) = 0.1; + fval = Inf; + exit_flag = 0; + return +end + +if imag(lik)~=0 + info(1) = 46; + info(2) = 0.1; + fval = Inf; + exit_flag = 0; + return +end + +% Add the (logged) prior density for the dsge-parameters. +lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4); +fval = (lik-lnprior); + +if isnan(fval) + info(1) = 47; + info(2) = 0.1; + fval = Inf; + exit_flag = 0; + return +end + +if imag(fval)~=0 + info(1) = 48; + info(2) = 0.1; + fval = Inf; + exit_flag = 0; + return +end + +if (nargout == 10) + if isinf(dsge_prior_weight) + iXX = iGXX; + else + iXX = tmp2; + end +end + +if (nargout==11) + 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*NumberOfObservations); + prior.DF = prior.ArtificialSampleSize - NumberOfParameters - NumberOfObservedVariables; + prior.iGXX = iGXX; +end + +if fval == Inf + pause +end \ No newline at end of file diff --git a/matlab/dsge_var_likelihood_1.m b/matlab/dsge_var_likelihood_1.m deleted file mode 100644 index e01764d26..000000000 --- a/matlab/dsge_var_likelihood_1.m +++ /dev/null @@ -1,312 +0,0 @@ -function [fval,info,exit_flag,grad,hess,SteadyState,trend_coeff,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood_1(xparam1,DynareDataset,DynareInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults) -% 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 SteadyState [double] Steady state vector possibly recomputed -% by call to dynare_results() -% o trend_coeff [double] place holder for trend coefficients, -% currently not supported by dsge_var -% 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. - -% Copyright (C) 2006-2012 Dynare Team -% -% 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 . - -persistent dsge_prior_weight_idx - -grad=[]; -hess=[]; -exit_flag = []; -info = 0; -PHI = []; -SIGMAu = []; -iXX = []; -prior = []; -SteadyState = []; -trend_coeff = []; - -% Initialization of of the index for parameter dsge_prior_weight in Model.params. -if isempty(dsge_prior_weight_idx) - dsge_prior_weight_idx = strmatch('dsge_prior_weight',Model.param_names); -end - -% Get the number of estimated (dsge) parameters. -nx = EstimatedParameters.nvx + EstimatedParameters.np; - -% Get the number of observed variables in the VAR model. -NumberOfObservedVariables = DynareDataset.vobs; - -% Get the number of observations. -NumberOfObservations = DynareDataset.nobs; - - -% Get the number of lags in the VAR model. -NumberOfLags = DynareOptions.dsge_varlag; - -% Get the number of parameters in the VAR model. -NumberOfParameters = NumberOfObservedVariables*NumberOfLags ; -if ~DynareOptions.noconstant - NumberOfParameters = NumberOfParameters + 1; -end - -% Get empirical second order moments for the observed variables. -mYY = evalin('base', 'mYY'); -mYX = evalin('base', 'mYX'); -mXY = evalin('base', 'mXY'); -mXX = evalin('base', 'mXX'); - -% Initialize some of the output arguments. -fval = []; -exit_flag = 1; - -% Return, with endogenous penalty, if some dsge-parameters are smaller than the lower bound of the prior domain. -if DynareOptions.mode_compute ~= 1 && any(xparam1 < BoundsInfo.lb) - fval = Inf; - exit_flag = 0; - info(1) = 41; - k = find(xparam1 < BoundsInfo.lb); - info(2) = sum((BoundsInfo.lb(k)-xparam1(k)).^2); - return; -end - -% Return, with endogenous penalty, if some dsge-parameters are greater than the upper bound of the prior domain. -if DynareOptions.mode_compute ~= 1 && any(xparam1 > BoundsInfo.ub) - fval = Inf; - exit_flag = 0; - info(1) = 42; - k = find(xparam1 > BoundsInfo.ub); - info(2) = sum((xparam1(k)-BoundsInfo.ub(k)).^2); - return; -end - -% Get the variance of each structural innovation. -Q = Model.Sigma_e; -for i=1:EstimatedParameters.nvx - k = EstimatedParameters.var_exo(i,1); - Q(k,k) = xparam1(i)*xparam1(i); -end -offset = EstimatedParameters.nvx; - -% Update Model.params and Model.Sigma_e. -Model.params(EstimatedParameters.param_vals(:,1)) = xparam1(offset+1:end); -Model.Sigma_e = Q; - -% Get the weight of the dsge prior. -dsge_prior_weight = Model.params(dsge_prior_weight_idx); - -% Is the dsge prior proper? -if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/ ... - NumberOfObservations; - fval = Inf; - exit_flag = 0; - info(1) = 51; - info(2) = abs(NumberOfObservations*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables)); - % info(2)=dsge_prior_weight; - % info(3)=(NumberOfParameters+NumberOfObservedVariables)/DynareDataset.nobs; - return -end - -%------------------------------------------------------------------------------ -% 2. call model setup & reduction program -%------------------------------------------------------------------------------ - -% Solve the Dsge model and get the matrices of the reduced form solution. T and R are the matrices of the -% state equation -[T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict'); - -% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol). -if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) == 8 || ... - info(1) == 22 || info(1) == 24 || info(1) == 25 || info(1) == 10 - fval = Inf; - info(2) = 0.1; - exit_flag = 0; - return -elseif info(1) == 3 || info(1) == 4 || info(1) == 19 || info(1) == 20 || info(1) == 21 - fval = Inf; - info(2) = 0.1; - exit_flag = 0; - return -end - -% Define the mean/steady state vector. -if ~DynareOptions.noconstant - if DynareOptions.loglinear - constant = transpose(log(SteadyState(BayesInfo.mfys))); - else - constant = transpose(SteadyState(BayesInfo.mfys)); - end -else - constant = zeros(1,NumberOfObservedVariables); -end - - -%------------------------------------------------------------------------------ -% 3. theoretical moments (second order) -%------------------------------------------------------------------------------ - -% Compute the theoretical second order moments -tmp0 = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); -mf = BayesInfo.mf1; - -% Get the non centered second order moments -TheoreticalAutoCovarianceOfTheObservedVariables = zeros(NumberOfObservedVariables,NumberOfObservedVariables,NumberOfLags+1); -TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) = tmp0(mf,mf)+constant'*constant; -for lag = 1:NumberOfLags - tmp0 = T*tmp0; - TheoreticalAutoCovarianceOfTheObservedVariables(:,:,lag+1) = tmp0(mf,mf) + constant'*constant; -end - -% Build the theoretical "covariance" between Y and X -GYX = zeros(NumberOfObservedVariables,NumberOfParameters); -for i=1:NumberOfLags - GYX(:,(i-1)*NumberOfObservedVariables+1:i*NumberOfObservedVariables) = TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1); -end -if ~DynareOptions.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 ~DynareOptions.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)% Evaluation of the likelihood of the dsge-var model when the dsge prior weight is finite. - tmp0 = dsge_prior_weight*NumberOfObservations*TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) + mYY ; - tmp1 = dsge_prior_weight*NumberOfObservations*GYX + mYX; - tmp2 = inv(dsge_prior_weight*NumberOfObservations*GXX+mXX); - SIGMAu = tmp0 - tmp1*tmp2*tmp1'; clear('tmp0'); - [SIGMAu_is_positive_definite, penalty] = ispd(SIGMAu); - if ~SIGMAu_is_positive_definite - fval = Inf; - info(1) = 52; - info(2) = penalty; - exit_flag = 0; - return; - end - SIGMAu = SIGMAu / (NumberOfObservations*(1+dsge_prior_weight)); - PHI = tmp2*tmp1'; clear('tmp1'); - prodlng1 = sum(gammaln(.5*((1+dsge_prior_weight)*NumberOfObservations- ... - NumberOfObservedVariables*NumberOfLags ... - +1-(1:NumberOfObservedVariables)'))); - prodlng2 = sum(gammaln(.5*(dsge_prior_weight*NumberOfObservations- ... - NumberOfObservedVariables*NumberOfLags ... - +1-(1:NumberOfObservedVariables)'))); - lik = .5*NumberOfObservedVariables*log(det(dsge_prior_weight*NumberOfObservations*GXX+mXX)) ... - + .5*((dsge_prior_weight+1)*NumberOfObservations-NumberOfParameters)*log(det((dsge_prior_weight+1)*NumberOfObservations*SIGMAu)) ... - - .5*NumberOfObservedVariables*log(det(dsge_prior_weight*NumberOfObservations*GXX)) ... - - .5*(dsge_prior_weight*NumberOfObservations-NumberOfParameters)*log(det(dsge_prior_weight*NumberOfObservations*(GYY-GYX*inv(GXX)*GYX'))) ... - + .5*NumberOfObservedVariables*NumberOfObservations*log(2*pi) ... - - .5*log(2)*NumberOfObservedVariables*((dsge_prior_weight+1)*NumberOfObservations-NumberOfParameters) ... - + .5*log(2)*NumberOfObservedVariables*(dsge_prior_weight*NumberOfObservations-NumberOfParameters) ... - - prodlng1 + prodlng2; -else% Evaluation of the likelihood of the dsge-var model when the dsge prior weight is infinite. - iGXX = inv(GXX); - SIGMAu = GYY - GYX*iGXX*transpose(GYX); - PHI = iGXX*transpose(GYX); - lik = NumberOfObservations * ( log(det(SIGMAu)) + NumberOfObservedVariables*log(2*pi) + ... - trace(inv(SIGMAu)*(mYY - transpose(mYX*PHI) - mYX*PHI + transpose(PHI)*mXX*PHI)/NumberOfObservations)); - lik = .5*lik;% Minus likelihood -end - -if isnan(lik) - info(1) = 45; - info(2) = 0.1; - fval = Inf; - exit_flag = 0; - return -end - -if imag(lik)~=0 - info(1) = 46; - info(2) = 0.1; - fval = Inf; - exit_flag = 0; - return -end - -% Add the (logged) prior density for the dsge-parameters. -lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4); -fval = (lik-lnprior); - -if isnan(fval) - info(1) = 47; - info(2) = 0.1; - fval = Inf; - exit_flag = 0; - return -end - -if imag(fval)~=0 - info(1) = 48; - info(2) = 0.1; - fval = Inf; - exit_flag = 0; - return -end - -if (nargout == 10) - if isinf(dsge_prior_weight) - iXX = iGXX; - else - iXX = tmp2; - end -end - -if (nargout==11) - 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*NumberOfObservations); - prior.DF = prior.ArtificialSampleSize - NumberOfParameters - NumberOfObservedVariables; - prior.iGXX = iGXX; -end - -if fval == Inf - pause -end \ No newline at end of file diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 9e3b89114..fe57db9c4 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -85,10 +85,10 @@ if ~options_.dsge_var error(['Estimation: Unknown filter ' options_.particle.filter_algorithm]) end else - objective_function = str2func('dsge_likelihood_1'); + objective_function = str2func('dsge_likelihood'); end else - objective_function = str2func('dsge_var_likelihood_1'); + objective_function = str2func('dsge_var_likelihood'); end [dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_, bayestopt_, bounds] = ... @@ -253,7 +253,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation if options_.analytic_derivation && strcmp(func2str(objective_function),'dsge_likelihood') options = options_.analytic_derivation; options.analytic_derivation = 2; - [junk1, junk2, junk3, junk4, hh] = feval(objective_function,xparam1, ... + [junk1, junk2, hh] = feval(objective_function,xparam1, ... dataset_,dataset_info,options_,M_, ... estim_params_,bayestopt_,bounds,oo_); elseif isequal(options_.mode_compute,4) || ... diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m index fe26bd35c..7eac2d199 100644 --- a/matlab/identification_analysis.m +++ b/matlab/identification_analysis.m @@ -148,7 +148,7 @@ if info(1)==0, dataset_ = dseries(oo_.endo_simul(options_.varobs_id,100+1:end)',dates('1Q1'), options_.varobs); derivatives_info.no_DLIK=1; bounds = prior_bounds(bayestopt_,options_); - [fval,info,cost_flag,DLIK,AHess,ys,trend_coeff,M_,options_,bayestopt_,oo_] = dsge_likelihood(params',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_,derivatives_info); + [fval,DLIK,AHess,cost_flag,ys,trend_coeff,info,M_,options_,bayestopt_,oo_] = dsge_likelihood(params',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_,derivatives_info); % fval = DsgeLikelihood(xparam1,data_info,options_,M_,estim_params_,bayestopt_,oo_); options_.analytic_derivation = analytic_derivation; AHess=-AHess; diff --git a/matlab/mode_check.m b/matlab/mode_check.m index 55a141438..4e37d41be 100644 --- a/matlab/mode_check.m +++ b/matlab/mode_check.m @@ -140,7 +140,7 @@ for plt = 1:nbplt, end for i=1:length(z) xx(kk) = z(i); - [fval, info, exit_flag] = feval(fun,xx,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults); + [fval, junk1, junk2, exit_flag] = feval(fun,xx,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults); if exit_flag y(i,1) = fval; else diff --git a/matlab/optimization/mr_hessian.m b/matlab/optimization/mr_hessian.m index c88e94cd7..0fab774ad 100644 --- a/matlab/optimization/mr_hessian.m +++ b/matlab/optimization/mr_hessian.m @@ -61,7 +61,7 @@ if init return end -[f0, exit_flag, ff0]=penalty_objective_function(x,func,penalty,varargin{:}); +[f0, ff0]=penalty_objective_function(x,func,penalty,varargin{:}); h2=varargin{7}.ub-varargin{7}.lb; hmax=varargin{7}.ub-x; hmax=min(hmax,x-varargin{7}.lb); @@ -93,7 +93,7 @@ while igtol && check==0 && jitgtol && check==0 && jit. -[loss,info,exit_flag,vx,junk]=osr_obj_1(x,i_params,i_var,weights); \ No newline at end of file +global M_ oo_ options_ optimal_Q_ it_ +% global ys_ Sigma_e_ endo_nbr exo_nbr optimal_Q_ it_ ykmin_ options_ + +junk = []; +exit_flag = 1; +vx = []; +info=0; +loss=[]; +% set parameters of the policiy rule +M_.params(i_params) = x; + +% don't change below until the part where the loss function is computed +it_ = M_.maximum_lag+1; +[dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_); + +switch info(1) + case 1 + loss = 1e8; + return + case 2 + loss = 1e8*min(1e3,info(2)); + return + case 3 + loss = 1e8*min(1e3,info(2)); + return + case 4 + loss = 1e8*min(1e3,info(2)); + return + case 5 + loss = 1e8; + return + case 6 + loss = 1e8*min(1e3,info(2)); + return + case 7 + loss = 1e8*min(1e3); + return + case 8 + loss = 1e8*min(1e3,info(2)); + return + case 9 + loss = 1e8*min(1e3,info(2)); + return + case 20 + loss = 1e8*min(1e3,info(2)); + return + otherwise + if info(1)~=0 + loss = 1e8; + return; + end +end + +vx = get_variance_of_endogenous_variables(dr,i_var); +loss = full(weights(:)'*vx(:)); diff --git a/matlab/osr_obj_1.m b/matlab/osr_obj_1.m deleted file mode 100644 index d889d3cee..000000000 --- a/matlab/osr_obj_1.m +++ /dev/null @@ -1,90 +0,0 @@ -function [loss,info,exit_flag,vx,junk]=osr_obj_1(x,i_params,i_var,weights) -% objective function for optimal simple rules (OSR) -% INPUTS -% x vector values of the parameters -% over which to optimize -% i_params vector index of optimizing parameters in M_.params -% i_var vector variables indices -% weights vector weights in the OSRs -% -% OUTPUTS -% loss scalar loss function returned to solver -% info vector info vector returned by resol -% exit_flag scalar exit flag returned to solver -% vx vector variances of the endogenous variables -% junk empty place holder for penalty_objective_function -% -% SPECIAL REQUIREMENTS -% none -% Copyright (C) 2005-2013 Dynare Team -% -% 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 . - -global M_ oo_ options_ optimal_Q_ it_ -% global ys_ Sigma_e_ endo_nbr exo_nbr optimal_Q_ it_ ykmin_ options_ - -junk = []; -exit_flag = 1; -vx = []; -info=0; -loss=[]; -% set parameters of the policiy rule -M_.params(i_params) = x; - -% don't change below until the part where the loss function is computed -it_ = M_.maximum_lag+1; -[dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_); - -switch info(1) - case 1 - loss = 1e8; - return - case 2 - loss = 1e8*min(1e3,info(2)); - return - case 3 - loss = 1e8*min(1e3,info(2)); - return - case 4 - loss = 1e8*min(1e3,info(2)); - return - case 5 - loss = 1e8; - return - case 6 - loss = 1e8*min(1e3,info(2)); - return - case 7 - loss = 1e8*min(1e3); - return - case 8 - loss = 1e8*min(1e3,info(2)); - return - case 9 - loss = 1e8*min(1e3,info(2)); - return - case 20 - loss = 1e8*min(1e3,info(2)); - return - otherwise - if info(1)~=0 - loss = 1e8; - return; - end -end - -vx = get_variance_of_endogenous_variables(dr,i_var); -loss = full(weights(:)'*vx(:)); diff --git a/tests/estimation/TaRB/fs2000_tarb.mod b/tests/estimation/TaRB/fs2000_tarb.mod index 2dc573d96..6586df195 100644 --- a/tests/estimation/TaRB/fs2000_tarb.mod +++ b/tests/estimation/TaRB/fs2000_tarb.mod @@ -120,4 +120,3 @@ tarb_mode_compute=4, tarb_new_block_probability=0.3, silent_optimizer ); -