Add routine for conditional likelihood (first order).
parent
644e72c33c
commit
b7693c3273
GPG Key ID:
295C1FE89E17EB3C
|
|
@ -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;
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
% 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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -1 +1 @@
|
|||
Subproject commit 0bc1539f45890f7df9ae70973396fa27b79cd293
|
||||
Subproject commit 4348f4d57f83aa8f96c3ec23899c33971e847428
|
||||
|
|
@ -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)))
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
Loading…
Reference in New Issue