From b7693c3273649b2d688b0f5233ff7ef560769067 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Ry=C3=BBk=29?= Date: Fri, 6 Jan 2023 11:54:24 +0100 Subject: [PATCH] Add routine for conditional likelihood (first order). --- matlab/default_option_values.m | 3 +- matlab/dsge_conditional_likelihood_1.m | 204 ++++++++++++++++++ matlab/dynare_estimation_1.m | 2 + preprocessor | 2 +- tests/Makefile.am | 4 + .../1/fs2000_estimation_conditional.mod | 92 ++++++++ .../1/fs2000_estimation_exact.mod | 85 ++++++++ 7 files changed, 390 insertions(+), 2 deletions(-) create mode 100644 matlab/dsge_conditional_likelihood_1.m create mode 100644 tests/estimation/conditional-likelihood/1/fs2000_estimation_conditional.mod create mode 100644 tests/estimation/conditional-likelihood/1/fs2000_estimation_exact.mod diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m index d44766494..de0114c59 100644 --- a/matlab/default_option_values.m +++ b/matlab/default_option_values.m @@ -780,4 +780,5 @@ options_.varobs_id=[]; %initialize field options_.pac.estimation.ols.share_of_optimizing_agents.lb = 0.0; options_.pac.estimation.ols.share_of_optimizing_agents.ub = 1.0; -end +options_.conditional_likelihood.status = false; +options_.conditional_likelihood.order = 1; \ No newline at end of file diff --git a/matlab/dsge_conditional_likelihood_1.m b/matlab/dsge_conditional_likelihood_1.m new file mode 100644 index 000000000..f8a16e4b8 --- /dev/null +++ b/matlab/dsge_conditional_likelihood_1.m @@ -0,0 +1,204 @@ +function [fval, info, exitflag, DLIK, Hess, SteadyState, trend_coeff, Model, DynareOptions, BayesInfo, DynareResults] = ... + dsge_conditional_likelihood_1(xparam1, DynareDataset, DatasetInfo, DynareOptions, Model, EstimatedParameters, BayesInfo, BoundsInfo, DynareResults, derivatives_info) + +% Copyright (C) 2017-2023 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 . + +% Initialization of the returned variables and others... +fval = []; +SteadyState = []; +trend_coeff = []; +exitflag = true; +info = zeros(4,1); +DLIK = []; +Hess = []; + +% Exit with error if analytical_derivation option is used. +if DynareOptions.analytic_derivation + error('The analytic_derivation and conditional_likelihood are not compatible!') +end + +% Ensure that xparam1 is a column vector. +% (Don't do the transformation if xparam1 is empty, otherwise it would become a +% 0×1 matrix, which create issues with older MATLABs when comparing with [] in +% check_bounds_and_definiteness_estimation) +if ~isempty(xparam1) + xparam1 = xparam1(:); +end + +%------------------------------------------------------------------------------ +% 1. Get the structural parameters & define penalties +%------------------------------------------------------------------------------ +Model = set_all_parameters(xparam1,EstimatedParameters,Model); + +[fval, info, exitflag, Q, H] = check_bounds_and_definiteness_estimation(xparam1, Model, EstimatedParameters, BoundsInfo); +if info(1) + return +end + +iQ_upper_chol = chol(inv(Q)); + +% Return an error if the interface for measurement errors is used. +if ~isequal(H, zeros(size(H))) || EstimatedParameters.ncn || EstimatedParameters.ncx + error('Option conditional_likelihood does not support declaration of measurement errors. You can specify the measurement errors in the model block directly by adding measurement equations.') +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, 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) + if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 ||... + info(1) == 20 || info(1) == 21 || info(1) == 23 || info(1) == 26 || ... + info(1) == 81 || info(1) == 84 || info(1) == 85 || info(1) == 86 || ... + info(1) == 401 || info(1) == 402 || info(1) == 403 || ... %cycle reduction + info(1) == 411 || info(1) == 412 || info(1) == 413 % logarithmic reduction + %meaningful second entry of output that can be used + fval = Inf; + info(4) = info(2); + exitflag = false; + return + else + fval = Inf; + info(4) = 0.1; + exitflag = false; + return + end +end + +% check endogenous prior restrictions +info = endogenous_prior_restrictions(T, R, Model, DynareOptions, DynareResults); +if info(1) + fval = Inf; + info(4)=info(2); + exitflag = false; + 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_addition, trend_coeff] = compute_trend_coefficients(Model, DynareOptions, DynareDataset.vobs, DynareDataset.nobs); + trend = repmat(constant, 1, DynareDataset.nobs)+trend_addition; +else + trend_coeff = zeros(DynareDataset.vobs, 1); + trend = repmat(constant, 1, DynareDataset.nobs); +end + +% Return an error if some observations are missing. +if DatasetInfo.missing.state + error('Option conditional_likelihood is not compatible with missing observations.') +end + +% Get the selection matrix (vector of row indices for T and R) +Z = BayesInfo.mf; + +% Get the number of observed variables. +pp = DynareDataset.vobs; + +% Get the number of variables in the state equations (state variables plus observed variables). +mm = size(T, 1); + +% Get the number of innovations. +rr = length(Q); + +% Return an error if the number of shocks is not equal to the number of observations. +if ~isequal(pp, rr) + error('With conditional_likelihood the number of innovations must be equal to the number of observed varilables!') +end + +% Remove the trend. +Y = transpose(DynareDataset.data)-trend; + +% Set state vector (deviation to steady state) +S = zeros(mm, 1); + +%------------------------------------------------------------------------------ +% 3. Evaluate the conditional likelihood +%------------------------------------------------------------------------------ + +Rtild = inv(R(Z,:)); +const = -.5*rr*log(2*pi); +const = const + log(abs(det(Rtild))) + sum(log(diag(iQ_upper_chol))); + +llik = zeros(size(Y, 2)); + +Ytild = Rtild*Y; +Ttild = Rtild*T(Z,:); + +for t=1:size(Y, 2) + epsilon = Ytild(:,t) - Ttild*S; + upsilon = iQ_upper_chol*epsilon; + S = T*S + R*epsilon; + llik(t) = const - .5*(upsilon'*upsilon); +end + +% Computes minus log-likelihood. +likelihood = -sum(llik(DynareOptions.presample+1:size(Y, 2))); + + +% ------------------------------------------------------------------------------ +% 5. Adds prior if necessary +% ------------------------------------------------------------------------------ + +lnprior = priordens(xparam1, BayesInfo.pshape, BayesInfo.p6, BayesInfo.p7, BayesInfo.p3, BayesInfo.p4); + +if DynareOptions.endogenous_prior==1 + [lnpriormom] = endogenous_prior(Y, Pstar, BayesInfo, H); + fval = (likelihood-lnprior-lnpriormom); +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(4) = 0.1; + exitflag = false; + return +end + +if imag(fval)~=0 + fval = Inf; + info(1) = 48; + info(4) = 0.1; + exitflag = false; + return +end \ No newline at end of file diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 9e8d1af5e..9de11f35b 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -98,6 +98,8 @@ if ~options_.dsge_var else if options_.occbin.likelihood.status && options_.occbin.likelihood.inversion_filter objective_function = str2func('occbin.IVF_posterior'); + elseif options_.conditional_likelihood.status && options_.conditional_likelihood.order==1 + objective_function = str2func('dsge_conditional_likelihood_1'); else objective_function = str2func('dsge_likelihood'); end diff --git a/preprocessor b/preprocessor index 0bc1539f4..4348f4d57 160000 --- a/preprocessor +++ b/preprocessor @@ -1 +1 @@ -Subproject commit 0bc1539f45890f7df9ae70973396fa27b79cd293 +Subproject commit 4348f4d57f83aa8f96c3ec23899c33971e847428 diff --git a/tests/Makefile.am b/tests/Makefile.am index 83b533806..478a2e4a4 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -81,6 +81,8 @@ MODFILES = \ estimation/method_of_moments/AFVRR/AFVRR_M0.mod \ estimation/method_of_moments/AFVRR/AFVRR_MFB.mod \ estimation/method_of_moments/AFVRR/AFVRR_MFB_RRA.mod \ + estimation/conditional-likelihood/1/fs2000_estimation_exact.mod \ + estimation/conditional-likelihood/1/fs2000_estimation_conditional.mod \ estimation/system_prior_restriction/Gali_2015.mod \ estimation/no_init_estimation_check_first_obs/fs2000_init_check.mod \ estimation/univariate/nls/staticmodel.mod \ @@ -1022,6 +1024,8 @@ solve_algo_12_14/purely_forward_12.o.trs: solve_algo_12_14/purely_forward_refere solve_algo_12_14/purely_forward_14.m.trs: solve_algo_12_14/purely_forward_reference.m.trs solve_algo_12_14/purely_forward_14.o.trs: solve_algo_12_14/purely_forward_reference.o.trs +estimation/conditional-likelihood/1/fs2000_estimation_conditional.m.trs: estimation/conditional-likelihood/1/fs2000_estimation_exact.m.trs +estimation/conditional-likelihood/1/fs2000_estimation_conditional.o.trs: estimation/conditional-likelihood/1/fs2000_estimation_exact.o.trs observation_trends_and_prefiltering/MCMC: m/observation_trends_and_prefiltering/MCMC o/observation_trends_and_prefiltering/MCMC m/observation_trends_and_prefiltering/MCMC: $(patsubst %.mod, %.m.trs, $(filter observation_trends_and_prefiltering/MCMC/%.mod, $(MODFILES))) diff --git a/tests/estimation/conditional-likelihood/1/fs2000_estimation_conditional.mod b/tests/estimation/conditional-likelihood/1/fs2000_estimation_conditional.mod new file mode 100644 index 000000000..724679ade --- /dev/null +++ b/tests/estimation/conditional-likelihood/1/fs2000_estimation_conditional.mod @@ -0,0 +1,92 @@ +var m P c e W R k d n l gy_obs gp_obs y dA; +varexo e_a e_m; + +parameters alp bet gam mst rho psi del; + +alp = 0.33; +bet = 0.99; +gam = 0.003; +mst = 1.011; +rho = 0.7; +psi = 0.787; +del = 0.02; + +model; +dA = exp(gam+e_a); +log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m; +-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0; +W = l/n; +-(psi/(1-psi))*(c*P/(1-n))+l/n = 0; +R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; +1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0; +c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1); +P*c = m; +m-1+d = l; +e = exp(e_a); +y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); +gy_obs = dA*y/y(-1); +gp_obs = (P/P(-1))*m(-1)/dA; +end; + +steady_state_model; + dA = exp(gam); + gst = 1/dA; + m = mst; + khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1)); + xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1); + nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp ); + n = xist/(nust+xist); + P = xist + nust; + k = khst*n; + + l = psi*mst*n/( (1-psi)*(1-n) ); + c = mst/P; + d = l - mst + 1; + y = k^alp*n^(1-alp)*gst^alp; + R = mst/bet; + W = l/n; + ist = y-c; + q = 1 - d; + + e = 1; + + gp_obs = m/dA; + gy_obs = dA; +end; + + +shocks; +var e_a; stderr 0.014; +var e_m; stderr 0.005; +end; + +steady; + +check; + +estimated_params; +alp, beta_pdf, 0.356, 0.02; +bet, beta_pdf, 0.993, 0.002; +gam, normal_pdf, 0.0085, 0.003; +mst, normal_pdf, 1.0002, 0.007; +rho, beta_pdf, 0.129, 0.223; +psi, beta_pdf, 0.65, 0.05; +del, beta_pdf, 0.01, 0.005; +stderr e_a, inv_gamma_pdf, 0.035449, inf; +stderr e_m, inv_gamma_pdf, 0.008862, inf; +end; + +varobs gp_obs gy_obs; + +options_.solve_tolf = 1e-12; + +tic +estimation(conditional_likelihood,order=1,datafile='../../fsdat_simul',nobs=192,mode_compute=4,loglinear,mh_replic=5000,mh_nblocks=2,mh_jscale=0.8); +toc + +exact_likelihood = load('fs2000_estimation_exact/Output/fs2000_estimation_exact_results.mat'); +oo_exact_likelihood = exact_likelihood.oo_; +oo_conditional_likelihood = oo_; + +diff = 100*(oo_conditional_likelihood.posterior.optimization.mode-oo_exact_likelihood.posterior.optimization.mode)./oo_exact_likelihood.posterior.optimization.mode; +diff diff --git a/tests/estimation/conditional-likelihood/1/fs2000_estimation_exact.mod b/tests/estimation/conditional-likelihood/1/fs2000_estimation_exact.mod new file mode 100644 index 000000000..33b230f52 --- /dev/null +++ b/tests/estimation/conditional-likelihood/1/fs2000_estimation_exact.mod @@ -0,0 +1,85 @@ +var m P c e W R k d n l gy_obs gp_obs y dA; +varexo e_a e_m; + +parameters alp bet gam mst rho psi del; + +alp = 0.33; +bet = 0.99; +gam = 0.003; +mst = 1.011; +rho = 0.7; +psi = 0.787; +del = 0.02; + +model; +dA = exp(gam+e_a); +log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m; +-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0; +W = l/n; +-(psi/(1-psi))*(c*P/(1-n))+l/n = 0; +R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; +1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0; +c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1); +P*c = m; +m-1+d = l; +e = exp(e_a); +y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); +gy_obs = dA*y/y(-1); +gp_obs = (P/P(-1))*m(-1)/dA; +end; + +steady_state_model; + dA = exp(gam); + gst = 1/dA; + m = mst; + khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1)); + xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1); + nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp ); + n = xist/(nust+xist); + P = xist + nust; + k = khst*n; + + l = psi*mst*n/( (1-psi)*(1-n) ); + c = mst/P; + d = l - mst + 1; + y = k^alp*n^(1-alp)*gst^alp; + R = mst/bet; + W = l/n; + ist = y-c; + q = 1 - d; + + e = 1; + + gp_obs = m/dA; + gy_obs = dA; +end; + + +shocks; +var e_a; stderr 0.014; +var e_m; stderr 0.005; +end; + +steady; + +check; + +estimated_params; +alp, beta_pdf, 0.356, 0.02; +bet, beta_pdf, 0.993, 0.002; +gam, normal_pdf, 0.0085, 0.003; +mst, normal_pdf, 1.0002, 0.007; +rho, beta_pdf, 0.129, 0.223; +psi, beta_pdf, 0.65, 0.05; +del, beta_pdf, 0.01, 0.005; +stderr e_a, inv_gamma_pdf, 0.035449, inf; +stderr e_m, inv_gamma_pdf, 0.008862, inf; +end; + +varobs gp_obs gy_obs; + +options_.solve_tolf = 1e-12; + +tic +estimation(order=1,datafile='../../fsdat_simul',nobs=192,mode_compute=4,loglinear,mh_replic=5000,mh_nblocks=2,mh_jscale=0.8); +toc