Good day to you all,

I am fairly new to “coding” and Dynare, so the following could be seen as a quite brutal, ugly and inefficient code. Basically, in the reproduction of a RBC model, I have to calibrate some parameters with respect to some empirical moments, namely,

```
% Target moment (to be matched)
eta_emp_moment = 0.0251 / 0.0136;
varphif_emp_moment = 0.4572 / 0.0136;
varphiB_emp_moment = 0.0150 / 0.0136;
kappa_emp_moment = 0.4058 / 0.0136;
```

They do not need to be these specifically but I show them for the sake of illustration. I am aware of other posts such as Perform moment matching with Matlab and Calibration to match second moments of real data, but I fail to see how to generalise them.

Allow me to present shortly the code, and you will understand my issue. My idea was to create a matrix with all the combinations of the parameters (in my case 4 parameters). In the fifth row of the matrix, I have zeros where the distance between the empirical and simulated moment will be stored. In a loop, I run Dynare feeding it all the combinations of the parameters, and computing the distance. In the end, I would find the combination of parameters with the shortest distance.

I do not have a close form solution of the steady state, only a sketch. Nonetheless, Dynare is able to run (when I impose some parameters according to some “educated” guesses). Concretely, the code is as follows:

## Loop moment matching idea

```
% Define ranges
eta_range = 1:0.1:3; % set of parameter values for eta
varphif_range = 0:0.1:0.5; % set of parameter values for varphif
varphiB_range = 0:0.1:0.5; % set of parameter values for varphiB
kappa_range = 0:0.1:0.5; % set of parameter values for kappa
% Initialize the result matrix
n = length(eta_range) * length(varphif_range) * length(varphiB_range) * length(kappa_range) ;
result_matrix = zeros(n, 4);
index = 1;
% Nested loops to generate the result matrix
for eta = eta_range
for varphif = varphif_range
for varphiB = varphiB_range
for kappa = kappa_range
% Assign the values to the result matrix
result_matrix(index, :) = [eta, varphif, varphiB, kappa];
index = index + 1;
end
end
end
end
result_matrix = result_matrix';
% We add a row of zeros where the check_dist will be stored.
result_matrix = vertcat(result_matrix, zeros(1, size(result_matrix, 2)));
for index = 1:length(result_matrix)
eta = result_matrix(1,index);
varphif = result_matrix(2,index);
varphiB = result_matrix(3,index);
kappa = result_matrix(4,index);
% ............ Steady state sketch ..............................
screen_S = screenbar; % arbitrary between 0 and 1
zBg_S = mu_Bg / (1 - rho_B);
zBb_S = mu_Bb / (1 - rho_B);
zK_S = 1;
Lambda_S = 1;
M_S = Gamma^(-1) * beta;
yk_S = Ykmean; % empricial average
yb_S = Ybmean; % empricial average
kk_S = Kkmean; % empricial average
kb_S = Kbmean; % empricial average
Y_S = yb_S + yk_S;
e_S = Kbmean * xi;
s_S = B + Kbmean - e_S;
A_S = s_S + e_S;
K_S = kk_S + kb_S;
rk_S = alpha * (yk_S / kk_S);
rgov_S = rbar;
r_S = rbar;
rb_S = v * yb_S / kb_S;
w_k_S = (1 - alpha) * (yk_S / H_f);
c_S = (1 / K_Y - Gamma + 1 - deltaT) * K_T;
pi_S = yb_S - deltaB * kb_S + rgov_S * B - r_S * s_S;
etilde_S = pi_S + e_S;
TR_S = omega1 * kb_S * exp(-omega2 + etilde_S / kb_S);
T_S = TR_S + (1 + rgov_S) * B - Gamma * B;
d_S = etilde_S + TR_S - Gamma * e_S - (phi / 2) * screen_S^2;
p_S = d_S * (M_S * Gamma / (1 - M_S * Gamma));
mu_S = Gamma - M_S * (1 + r_S) * (1 - omega1 * TR_S / kb_S);
n1 = c_S + Gamma * s_S + Gamma * kk_S + p_S;
n2 = d_S + p_S + (1 + r_S) * s_S + (w_k_S * H_f) + (rk_S + 1 - deltaK) * kk_S - T_S;
n_S = min([n1, n2]);
inv_k_S = kk_S - (1 - deltaK) * kk_S / Gamma;
inv_b_S = kb_S - (1 - deltaB) * kb_S / Gamma;
% Save in param_km all the parameters and steady states
save param_km Gamma rho_f sigma_f H_f B ...
rho_B sigma_fB alpha beta v ...
xi theta omega1 omega2 eta ...
varphif varphiB kappa mu_Bg mu_Bb ...
phi deltaB deltaK Effectb Effectk ...
mu_K sigmaBg sigmaBb dbar ...
Y_S yk_S yb_S A_S s_S ...
T_S e_S c_S d_S kk_S ...
kb_S K_S TR_S M_S n_S ...
etilde_S p_S r_S rk_S rb_S ...
w_k_S rgov_S mu_S screen_S pi_S ...
zBg_S zBb_S zK_S Lambda_S inv_k_S ...
inv_b_S;
% ............ Inovoke dynare ..............................
dynare RBCmodel.mod noclearall nolog nograph
% if info(1)~=0
% error('Simulation failed')
% end
% ............ Moments of interest ..............................
% Y
Y_simul = oo_.endo_simul(strmatch('Y',M_.endo_names,'exact'),:); % Position of Y
std_Y_simul= std(diff(Y_simul) ./ Y_simul(1:end-1));
% s_c
s_simul = oo_.endo_simul(strmatch('s',M_.endo_names,'exact'),:); % Position of s
c_simul = oo_.endo_simul(strmatch('c',M_.endo_names,'exact'),:); % Position of c
s_c_simul_vec = s_simul./c_simul;
std_s_c_simul_vec = std(diff(s_c_simul_vec) ./ s_c_simul_vec(1:end-1));
% kk
kk_simul = oo_.endo_simul(strmatch('kk',M_.endo_names,'exact'),:); % Position of kk
std_kk_simul_vec = std(diff(kk_simul) ./ kk_simul(1:end-1));
% kb
kb_simul = oo_.endo_simul(strmatch('kb',M_.endo_names,'exact'),:); % Position of kb
std_kb_simul_vec = std(diff(kb_simul) ./ kb_simul(1:end-1));
% d
d_simul = oo_.endo_simul(strmatch('d',M_.endo_names,'exact'),:); % Position of d
std_d_simul_vec = std(diff(d_simul) ./ d_simul(1:end-1));
% Target eta
eta_sim_moment = std_s_c_simul_vec/std_Y_simul;
% Target varphif
varphif_sim_moment = std_kk_simul_vec/std_Y_simul;
% Target varphiB
varphiB_sim_moment = std_kb_simul_vec/std_Y_simul;
% Target kappa
kappa_sim_moment = std_d_simul_vec/std_Y_simul;
% ............ Computation of check_dist ........................
if info(1)==0
eta_dist_moment = (eta_emp_moment - eta_sim_moment)^2 ;
varphif_dist_moment = (varphif_emp_moment - varphif_sim_moment)^2 ;
varphiB_dist_moment = (varphiB_emp_moment - varphiB_sim_moment)^2 ;
kappa_dist_moment = (kappa_emp_moment - kappa_sim_moment)^2 ;
% Compute distance
check_dist = eta_dist_moment + varphif_dist_moment + varphiB_dist_moment + kappa_dist_moment;
result_matrix(5,index) = check_dist ;
else
check_dist = 1000 ; % or NA (because combination has no result)
result_matrix(5,index) = check_dist ;
end
end
```

This loop does not run because I do not know how to tell dynare to flag when it was not able to estimate the steady state and jump to the next iteration. I read in other posts that `info(1)`

was used but with `stoch_simul(M_, options_, oo_, var_list_)`

and that looping over dynare is very inefficient, but I do not know if `stoch_simul()`

re-estimates the steady state. Is there a way of solving this problem?

I thank in advance for any help or comment, and sorry for the long and most likely confusing post.

Diego H.