Matching multiple parameters to target moments

Steady state
1

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.

2

Dear prof. @jpfeifer , can you confirm that, in order to loop over parameters (for matching moments, finding optimal welfare, etc.),

  1. an external steadystate m file needs to be created (as shown by willi mutschler on youtube and in other posts);
  2. the external steadystate m can give an approximation of the steady states (it does not need to be precise because dynare will estimate the steady states anyways);
  3. after having invoked dynare in the first loop, you need to use [info, oo_, options_, M_] = stoch_simul(M_, options_, oo_, var_list_) in the rest of the loops, as shown for example in Parameter loop in Dynare 5.x - #3 by jpfeifer;
  4. after that, you can retrieve the information needed using M_ and oo_

Sorry for the inconvenience and I thank you in advance!

3
  1. No, you don’t need an external steady state file, but it may be very useful to efficiently solve the problem.
  2. Steady state files need to be exact. Otherwise, we are just talking initial values.
  3. Yes, ideally that is the way to proceed as it will be faster. stoch_simul will recompute the steady state.
  4. Is there a reason why you do a grid-search type matching instead of e.g. doing GMM or SMM estimation in Dynare.
  5. I cannot run your code due to

Unrecognized function or variable ‘screenbar’.

4

I cannot thank you enough for your answers!

  1. I see! For my model, I do not have the exact steady state, just an approximation, as shown below % ... Steady state sketch ... in the code presented. Provided I have sound parameters, dynare is able to run.

  2. In my case, in param_km steady state values and parameters are saved. Then I feed them to the RBCmodel.mod file which estimates the steady states (inside I have load(param_km)). So I am basically taking initial values.

  3. Ok!

  4. I try to follow what the paper did… It is tricky because parameters and steady states are intertwined:

    • some parameters match data moments,
    • others are selected using steady-state conditions (I solve a system of equations composed by a subset of steady state equations which depend on data), and
    • the rest depend on simulated moments (my current case with the loop which selects a combination of parameters that reduces the distance between data and simulated moments).

    I wish I could use GMM, SMM or Bayesian estimation, but you need to have an analytical solution for the steady state, would you not?

  5. The code is just an outline. If needed, I can provide it.

To sum up: is it possible to loop over stoch_simul using only initial values for the steady states as in the code presented? I was trying to find a way of ignoring errors of the dynare command when it could not find steady states and jump to the next iteration.

5

You should aim for an analytical steady state. That does not mean computing it in closed form. Most of the time, you reduce the problem to a few equations and then solve them numerically.
Just discarding draws where the steady state cannot be computed is dangerous because you may erroneously discard draws that fit the data very well.

That being said, you can run stoch_simul using only initial values.

6

I know some steady state values (for example, Lambda_S = 1, rgov_S = r_S or rk_S = (Gamma/beta + 1 + deltaK), and others), but it is not enough. For capital, for example, I used kk_S = Kkmean; and kb_S = Kbmean as initial values (drawn from data).

  1. In the mod file, should I create a block with steady_state_model for the variables I know and initval block for the variables I do not know but have an approximation?
  2. I still do not understand the loop according to steady state guesses. I was looking at Parameter loop in Dynare 5.x - #3 by jpfeifer, Loop over parameters - #67 by cmarch and Looping over initial values - guide/example code - #6 by jpfeifer (with code matlab/evaluate_steady_state.m · master · Dynare / dynare · GitLab), to test the ideas presented. I tried to simplify it to just see the changes in the steady state value of consumption c, but c_storage remains empty. That loop starts in line 414 of Loop_ss_eval.m. In line 291, I have the sketch of steady states and, in 411, I invoke dynare.
    SteadyStateLoop.zip (45.7 KB)
7
  1. No, it’s either steady_state_model or initval. You cannot do both.
  2. Please use Dynare 6 and
        [steady_state,params,info] = evaluate_steady_state(oo_.steady_state,oo_.exo_steady_state,M_,options_,~options_.steadystate.nocheck);

The correctness check needs to be

if info(1)==0
8

Thank you! Just to confirm, when you mention the function evaluate_steady_state, are you referring to matlab/evaluate_steady_state.m · master · Dynare / dynare · GitLab (or is there another function because I looked at the manual and it is no where to be found)? If so,

  1. I noticed that from line 490 onwards (and in other places throughout the code), seemingly external functions are included. In the aforementioned webpage, should I find them and download them as well?

  2. does the evaluate_steady_state function also evaluates if the given parameters violate the Blanchard Kahn conditions?

9
  1. The evaluate_steady_state function is contained in every Dynare distribution and is the function Dynare internally uses to compute the steady state values. There is no reason to download it again. We do not put it in the manual, because usually users don’t need to call it manually.
  2. Yes, that function only checks the steady state. If you also need the Blanchard Kahn conditions, you would call resol.
10

I see! Thank you!
As a follow up question, I am looking at your resol function (matlab/stochastic_solver/resol.m · master · Dynare / dynare · GitLab). In it, one of the outputs is params and the description states: vector of potentially updated parameters. I then ask: can the function solve for a parameter? If so, how could I that? Via the function or would I need to code a loop or something?

11

resol will potentially call a steady state file, which in turn can updated parameters. That’s the reason for the phrasing. You would still need to pass all necessary parameters via M_.params..