Calibration / method of moments - Gali (2015)
This tutorial is intended to show the workflow to calibrate a model using the method of moments. The tutorial is based on a standard model of monetary policy and will showcase the the use of gradient based optimisers and 2nd and 3rd order pruned solutions.
Define the model
The first step is always to name the model and write down the equations. For the Galı́ (2015), Chapter 3 this would go as follows:
julia> using MacroModellingjulia> @model Gali_2015 begin W_real[0] = C[0] ^ σ * N[0] ^ φ Q[0] = β * (C[1] / C[0]) ^ (-σ) * Z[1] / Z[0] / Pi[1] R[0] = 1 / Q[0] Y[0] = A[0] * (N[0] / S[0]) ^ (1 - α) R[0] = Pi[1] * realinterest[0] R[0] = 1 / β * Pi[0] ^ ϕᵖⁱ * (Y[0] / Y[ss]) ^ ϕʸ * exp(nu[0]) C[0] = Y[0] log(A[0]) = ρ_a * log(A[-1]) + std_a * eps_a[x] log(Z[0]) = ρ_z * log(Z[-1]) - std_z * eps_z[x] nu[0] = ρ_ν * nu[-1] + std_nu * eps_nu[x] MC[0] = W_real[0] / (S[0] * Y[0] * (1 - α) / N[0]) 1 = θ * Pi[0] ^ (ϵ - 1) + (1 - θ) * Pi_star[0] ^ (1 - ϵ) S[0] = (1 - θ) * Pi_star[0] ^ (( - ϵ) / (1 - α)) + θ * Pi[0] ^ (ϵ / (1 - α)) * S[-1] Pi_star[0] ^ (1 + ϵ * α / (1 - α)) = ϵ * x_aux_1[0] / x_aux_2[0] * (1 - τ) / (ϵ - 1) x_aux_1[0] = MC[0] * Y[0] * Z[0] * C[0] ^ (-σ) + β * θ * Pi[1] ^ (ϵ + α * ϵ / (1 - α)) * x_aux_1[1] x_aux_2[0] = Y[0] * Z[0] * C[0] ^ (-σ) + β * θ * Pi[1] ^ (ϵ - 1) * x_aux_2[1] log_y[0] = log(Y[0]) log_W_real[0] = log(W_real[0]) log_N[0] = log(N[0]) pi_ann[0] = 4 * log(Pi[0]) i_ann[0] = 4 * log(R[0]) r_real_ann[0] = 4 * log(realinterest[0]) M_real[0] = Y[0] / R[0] ^ η endModel: Gali_2015 Variables Total: 23 Auxiliary: 0 States: 4 Auxiliary: 0 Jumpers: 5 Auxiliary: 0 Shocks: 3 Parameters: 16
First, we load the package and then use the @model macro to define our model. The first argument after @model is the model name and will be the name of the object in the global environment containing all information regarding the model. The second argument to the macro are the equations, which we write down between begin and end. Equations can contain an equality sign or the expression is assumed to equal 0. Equations cannot span multiple lines (unless you wrap the expression in brackets) and the timing of endogenous variables are expressed in the square brackets following the variable name (e.g. [-1] for the past period). Exogenous variables (shocks) are followed by a keyword in square brackets indicating them being exogenous (in this case [x]). Note that names can leverage julia's unicode capabilities (e.g. alpha can be written as α).
Define the parameters
Next we need to add the parameters of the model. The macro @parameters takes care of this:
julia> @parameters Gali_2015 begin σ = 1 φ = 5 ϕᵖⁱ = 1.5 ϕʸ = 0.125 θ = 0.75 ρ_ν = 0.5 ρ_z = 0.5 ρ_a = 0.9 β = 0.99 η = 3.77 α = 0.25 ϵ = 9 τ = 0 std_a = .01 std_z = .05 std_nu = .0025 endRemove redundant variables in non stochastic steady state problem: 1.394 seconds Set up non stochastic steady state problem: 0.592 seconds Take symbolic derivatives up to first order: 0.407 seconds Find non stochastic steady state: 6.449 seconds Model: Gali_2015 Variables Total: 23 Auxiliary: 0 States: 4 Auxiliary: 0 Jumpers: 5 Auxiliary: 0 Shocks: 3 Parameters: 16
The block defining the parameters above only describes the simple parameter definitions the same way you assign values (e.g. α = .25).
Note that we have to write one parameter definition per line.
Linear solution
Inspect model moments
Given the equations and parameters, we have everything to we need for the package to generate the theoretical model moments. You can retrieve the mean of the linearised model as follows:
julia> get_mean(Gali_2015)1-dimensional KeyedArray(NamedDimsArray(...)) with keys: ↓ Variables ∈ 23-element Vector{Symbol} And data, 23-element Vector{Float64}: (:A) 1.0 (:C) 0.9505798249541407 (:MC) 0.8888888888888887 (:M_real) 0.9152363832868945 (:N) 0.934655265184067 (:Pi) 0.9999999999999991 (:Pi_star) 0.9999999999999973 (:Q) 0.9900000000000011 ⋮ (:log_y) -0.05068313851352055 (:nu) 0.0 (:pi_ann) -3.5527136788005025e-15 (:r_real_ann) 0.04020134341400426 (:realinterest) 1.0101010101010097 (:x_aux_1) 3.451995685005287 (:x_aux_2) 3.8834951456309885
and the standard deviation like this:
julia> get_standard_deviation(Gali_2015)1-dimensional KeyedArray(NamedDimsArray(...)) with keys: ↓ Variables ∈ 23-element Vector{Symbol} And data, 23-element Vector{Float64}: (:A) 0.022941573387055825 (:C) 0.03357169677998865 (:MC) 0.21609085046272047 (:M_real) 0.059266198320989456 (:N) 0.03786945958932778 (:Pi) 0.012358762176559656 (:Pi_star) 0.037076286529677975 (:Q) 0.02046853220935804 ⋮ (:log_y) 0.035317072694666385 (:nu) 0.0028867513459481255 (:pi_ann) 0.04943504870623866 (:r_real_ann) 0.05644645066322982 (:realinterest) 0.01425415420788631 (:x_aux_1) 0.9515263512142034 (:x_aux_2) 0.5166586569325933
You could also simply use: std or get_std to the same effect.
Another interesting output is the autocorrelation of the model variables:
julia> get_autocorrelation(Gali_2015)2-dimensional KeyedArray(NamedDimsArray(...)) with keys: ↓ Variables ∈ 23-element Vector{Symbol} → Autocorrelation_orders ∈ 5-element UnitRange{Int64} And data, 23×5 Matrix{Float64}: (1) (2) … (4) (5) (:A) 0.9 0.81 0.6561 0.59049 (:C) 0.610108 0.404152 0.225901 0.185193 (:MC) 0.508432 0.261805 0.0750132 0.0430389 (:M_real) 0.729895 0.571853 0.403664 0.352666 (:N) 0.508432 0.261805 … 0.0750132 0.0430389 (:Pi) 0.626445 0.427023 0.250145 0.208033 ⋮ ⋱ ⋮ (:log_y) 0.610108 0.404152 0.225901 0.185193 (:nu) 0.5 0.25 0.0625 0.03125 (:pi_ann) 0.626445 0.427023 … 0.250145 0.208033 (:r_real_ann) 0.506897 0.259655 0.0727346 0.0408922 (:realinterest) 0.506897 0.259655 0.0727346 0.0408922 (:x_aux_1) 0.700916 0.531282 0.360659 0.31215 (:x_aux_2) 0.783352 0.646693 0.482995 0.427405
or the covariance:
julia> get_covariance(Gali_2015)2-dimensional KeyedArray(NamedDimsArray(...)) with keys: ↓ Variables ∈ 23-element Vector{Symbol} → 𝑉𝑎𝑟𝑖𝑎𝑏𝑙𝑒𝑠 ∈ 23-element Vector{Symbol} And data, 23×23 Matrix{Float64}: (:A) (:C) … (:x_aux_1) (:x_aux_2) (:A) 0.000526316 0.000404089 -0.0154711 -0.00997609 (:C) 0.000404089 0.00112706 0.00730576 0.000310327 (:MC) -0.000719776 0.0055578 0.16467 0.0732626 (:M_real) 0.00103078 -0.000276243 -0.0554553 -0.0300455 (:N) -0.000126139 0.000973992 … 0.028858 0.0128391 (:Pi) -0.000159411 0.000169707 0.0115462 0.00587158 ⋮ ⋱ (:log_y) 0.000425097 0.00118565 0.00768558 0.000326461 (:nu) -1.68584e-17 -8.20937e-6 -0.00016948 -5.38539e-5 (:pi_ann) -0.000637646 0.000678826 … 0.0461849 0.0234863 (:r_real_ann) -0.000170039 0.00143464 0.0418524 0.0185996 (:realinterest) -4.29391e-5 0.000362283 0.0105688 0.00469687 (:x_aux_1) -0.0154711 0.00730576 0.905402 0.480501 (:x_aux_2) -0.00997609 0.000310327 0.480501 0.266936
Parameter sensitivities
Before embarking on calibrating the model it is useful to get familiar with the impact of parameter changes on model moments. MacroModelling.jl provides the partial derivatives of the model moments with respect to the model parameters. The model we are working with is of a medium size and by default derivatives are automatically shown as long as the calculation does not take too long (too many derivatives need to be taken). In this case they are not shown but it is possible to show them by explicitly defining the parameter for which to take the partial derivatives for:
julia> get_mean(Gali_2015, parameter_derivatives = :σ)2-dimensional KeyedArray(NamedDimsArray(...)) with keys: ↓ Variables ∈ 23-element Vector{Symbol} → Mean_and_∂mean∂parameter ∈ 2-element Vector{Symbol} And data, 23×2 Matrix{Float64}: (:Mean) (:σ) (:A) 1.0 0.0 (:C) 0.95058 0.0060223 (:MC) 0.888889 -9.54098e-18 (:M_real) 0.915236 0.00579838 (:N) 0.934655 0.00789521 (:Pi) 1.0 -0.0 ⋮ (:log_y) -0.0506831 0.00633539 (:nu) 0.0 0.0 (:pi_ann) -3.55271e-15 -0.0 (:r_real_ann) 0.0402013 0.0 (:realinterest) 1.0101 0.0 (:x_aux_1) 3.452 0.174958 (:x_aux_2) 3.8835 0.196828
or for multiple parameters:
julia> get_mean(Gali_2015, parameter_derivatives = [:σ, :α, :β, :ϕᵖⁱ, :φ])2-dimensional KeyedArray(NamedDimsArray(...)) with keys: ↓ Variables ∈ 23-element Vector{Symbol} → Mean_and_∂mean∂parameter ∈ 6-element Vector{Symbol} And data, 23×6 Matrix{Float64}: (:Mean) (:σ) … (:β) (:α) (:A) 1.0 0.0 0.0 0.0 (:C) 0.95058 0.0060223 4.78473e-15 -0.0941921 (:MC) 0.888889 -9.54098e-18 3.60294e-14 -1.33227e-15 (:M_real) 0.915236 0.00579838 3.48529 -0.09069 (:N) 0.934655 0.00789521 … 6.3651e-15 -0.207701 (:Pi) 1.0 -0.0 -1.46549e-16 -0.0 ⋮ ⋱ ⋮ (:log_y) -0.0506831 0.00633539 5.03348e-15 -0.0990891 (:nu) 0.0 0.0 0.0 0.0 (:pi_ann) -3.55271e-15 -0.0 … -5.86198e-16 -0.0 (:r_real_ann) 0.0402013 0.0 -4.0404 0.0 (:realinterest) 1.0101 0.0 -1.0203 0.0 (:x_aux_1) 3.452 0.174958 10.0544 -1.47167e-13 (:x_aux_2) 3.8835 0.196828 11.3112 -0.0
We can do the same for standard deviation or variance, and all parameters:
julia> get_std(Gali_2015, parameter_derivatives = get_parameters(Gali_2015))2-dimensional KeyedArray(NamedDimsArray(...)) with keys: ↓ Variables ∈ 23-element Vector{Symbol} → Standard_deviation_and_∂standard_deviation∂parameter ∈ 17-element Vector{Symbol} And data, 23×17 Matrix{Float64}: (:Standard_deviation) … (:std_z) (:std_nu) (:A) 0.0229416 -9.40441e-34 5.48591e-34 (:C) 0.0335717 0.48179 0.0963579 (:MC) 0.216091 4.18882 0.837765 (:M_real) 0.0592662 0.491332 0.254845 (:N) 0.0378695 … 0.734082 0.146816 (:Pi) 0.0123588 0.167366 0.0334732 ⋮ ⋱ (:log_y) 0.0353171 0.506838 0.101368 (:nu) 0.00288675 0.0 1.1547 (:pi_ann) 0.049435 … 0.669465 0.133893 (:r_real_ann) 0.0564465 1.09678 0.253692 (:realinterest) 0.0142542 0.276965 0.0640637 (:x_aux_1) 0.951526 9.39926 1.44896 (:x_aux_2) 0.516659 2.99988 0.269446
julia> get_variance(Gali_2015, parameter_derivatives = get_parameters(Gali_2015))2-dimensional KeyedArray(NamedDimsArray(...)) with keys: ↓ Variables ∈ 23-element Vector{Symbol} → Variance_and_∂variance∂parameter ∈ 17-element Vector{Symbol} And data, 23×17 Matrix{Real}: (:Variance) (:σ) … (:std_z) (:std_nu) (:A) 0.000526316 -9.54188e-20 -4.31504e-35 2.51711e-35 (:C) 0.00112706 -0.00101983 0.032349 0.0064698 (:MC) 0.0466953 -0.0394622 1.81033 0.362067 (:M_real) 0.00351248 -0.000388635 0.0582387 0.0302074 (:N) 0.0014341 -0.00150695 … 0.0555986 0.0111197 (:Pi) 0.000152739 -6.64228e-5 0.00413688 0.000827376 ⋮ ⋱ (:log_y) 0.0012473 -0.00114444 0.0358 0.00716001 (:nu) 8.33333e-6 -6.68083e-22 0.0 0.00666667 (:pi_ann) 0.00244382 -0.00106276 … 0.06619 0.013238 (:r_real_ann) 0.0031862 -0.00279404 0.123819 0.02864 (:realinterest) 0.000203181 -0.000178173 0.00789579 0.00182635 (:x_aux_1) 0.905402 -0.00995817 17.8873 2.75744 (:x_aux_2) 0.266936 0.100826 3.09983 0.278423
You can use this information to calibrate certain values to your targets. For example, let's say we want to have higher real wages (:W_real), and lower inflation volatility. Since there are too many variables and parameters for them to be shown here, let's print only a subset of them:
julia> get_mean(Gali_2015, parameter_derivatives = [:σ, :std_a, :α], variables = [:W_real,:Pi])2-dimensional KeyedArray(NamedDimsArray(...)) with keys: ↓ Variables ∈ 2-element Vector{Symbol} → Mean_and_∂mean∂parameter ∈ 4-element Vector{Symbol} And data, 2×4 Matrix{Float64}: (:Mean) (:σ) (:α) (:std_a) (:Pi) 1.0 -0.0 -0.0 -0.0 (:W_real) 0.678025 -0.00143185 -0.820546 -0.0
julia> get_std(Gali_2015, parameter_derivatives = [:σ, :std_a, :α], variables = [:W_real,:Pi])2-dimensional KeyedArray(NamedDimsArray(...)) with keys: ↓ Variables ∈ 2-element Vector{Symbol} → Standard_deviation_and_∂standard_deviation∂parameter ∈ 4-element Vector{Symbol} And data, 2×4 Matrix{Float64}: (:Standard_deviation) (:σ) (:α) (:std_a) (:Pi) 0.0123588 -0.00268728 -0.0166021 0.390677 (:W_real) 0.156462 -0.0674815 0.141894 0.0348056
Looking at the sensitivity table we see that lowering the production function parameter :α will increase real wages, but at the same time it will increase inflation volatility. We could compensate that effect by decreasing the standard deviation of the total factor productivity shock :std_a.
Method of moments
Instead of doing this by hand we can also set a target and have an optimiser find the corresponding parameter values. In order to do that we need to define targets, and set up an optimisation problem.
Our targets are:
- Mean of
W_real = 0.7 - Standard deviation of
Pi = 0.01
For the optimisation problem we use the L-BFGS algorithm implemented in Optim.jl. This optimisation algorithm is very efficient and gradient based. Note that all model outputs are differentiable with respect to the parameters using automatic and implicit differentiation.
The package provides functions specialised for the use with gradient based code (e.g. gradient-based optimisers or samplers). For model statistics we can use get_statistics to get the mean of real wages and the standard deviation of inflation like this:
julia> get_statistics(Gali_2015, Gali_2015.parameter_values, parameters = Gali_2015.parameters, mean = [:W_real], standard_deviation = [:Pi])2-element Vector{AbstractArray{Float64}}: [0.6780252644037245] [0.012358762176559656]
First we pass on the model object, followed by the parameter values and the parameter names the values correspond to. Then we define the outputs we want: for the mean we want real wages and for the standard deviation we want inflation. We can also get outputs for variance, covariance, or autocorrelation the same way as for the mean and standard deviation.
Next, let's define a function measuring how close we are to our target for given values of :α and :std_a:
julia> function distance_to_target(parameter_value_inputs) model_statistics = get_statistics(Gali_2015, parameter_value_inputs, parameters = [:α, :std_a], mean = [:W_real], standard_deviation = [:Pi]) targets = [0.7, 0.01] return sum(abs2, vcat(model_statistics...) - targets) enddistance_to_target (generic function with 1 method)
Now let's test the function with the current parameter values. In case we forgot the parameter values we can also look them up like this:
julia> get_parameters(Gali_2015, values = true)16-element Vector{Pair{String, Float64}}: "σ" => 1.0 "φ" => 5.0 "ϕᵖⁱ" => 1.5 "ϕʸ" => 0.125 "θ" => 0.75 "ρ_ν" => 0.5 "ρ_z" => 0.5 "ρ_a" => 0.9 "β" => 0.99 "η" => 3.77 "α" => 0.25 "ϵ" => 9.0 "τ" => 0.0 "std_a" => 0.01 "std_z" => 0.05 "std_nu" => 0.0025
with this we can test the distance function:
julia> distance_to_target([0.25, 0.01])0.0004884527635317859
Next we can pass it on to an optimiser and find the parameters corresponding to the best fit like this:
julia> using Optim, LineSearchesjulia> sol = Optim.optimize(distance_to_target, [0,0], [1,1], [0.25, 0.01], Optim.Fminbox(Optim.LBFGS(linesearch = LineSearches.BackTracking(order = 3))))* Status: success * Candidate solution Final objective value: 5.107257e-07 * Found with Algorithm: Fminbox with L-BFGS * Convergence measures |x - x'| = 2.29e-06 ≰ 0.0e+00 |x - x'|/|x'| = 1.03e-05 ≰ 0.0e+00 |f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00 |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00 |g(x)| = 8.72e-09 ≤ 1.0e-08 * Work counters Seconds run: 1 (vs limit Inf) Iterations: 4 f(x) calls: 26 ∇f(x) calls: 24
The first argument to the optimisation call is the function we defined previously, followed by lower and upper bounds, the starting values, and finally the algorithm. For the algorithm we have to add Fminbox because we have bounds (optional) and we set the specific line search method to speed up convergence (recommended but optional).
The output shows that we could almost perfectly match the target and the values of the parameters found by the optimiser are:
julia> sol.minimizer2-element Vector{Float64}: 0.223302559538744 9.506455316606803e-8
slightly lower for both parameters (in line with what we understood from the sensitivities).
You can combine the method of moments with estimation by simply adding the distance to the target to the posterior loglikelihood.
Nonlinear solutions
So far we used the linearised solution of the model. The package also provides nonlinear solutions and can calculate the theoretical model moments for pruned second and third order perturbation solutions. This can be of interest because nonlinear solutions capture volatility effects (at second order) and asymmetries (at third order). Furthermore, the moments of the data are often non-gaussian while linear solutions with gaussian noise can only generate gaussian distributions of model variables. Nonetheless, already pruned second order solutions produce non-gaussian skewness and kurtosis with gaussian noise.
From a user perspective little changes other than specifying that the solution algorithm is :pruned_second_order or :pruned_third_order.
For example we can get the mean for the pruned second order solution:
julia> get_mean(Gali_2015, parameter_derivatives = [:σ, :std_a, :α], variables = [:W_real,:Pi], algorithm = :pruned_second_order)Take symbolic derivatives up to second order: 1.127 seconds 2-dimensional KeyedArray(NamedDimsArray(...)) with keys: ↓ Variables ∈ 2-element Vector{Symbol} → Mean_and_∂mean∂parameter ∈ 4-element Vector{Symbol} And data, 2×4 Matrix{Float64}: (:Mean) (:σ) (:α) (:std_a) (:Pi) 1.00915 -0.00566878 0.0136518 0.509002 (:W_real) 0.650863 0.0135465 -0.823786 -1.15286
Note that the mean of real wages is lower, while inflation is higher. We can see the effect of volatility with the partial derivatives for the shock standard deviations being non-zero. Larger shocks sizes drive down the mean of real wages while they increase inflation.
The mean of the variables does not change if we use pruned third order perturbation by construction but the standard deviation does. Let's look at the standard deviations for the pruned second order solution first:
julia> get_std(Gali_2015, parameter_derivatives = [:σ, :std_a, :α], variables = [:W_real,:Pi], algorithm = :pruned_second_order)2-dimensional KeyedArray(NamedDimsArray(...)) with keys: ↓ Variables ∈ 2-element Vector{Symbol} → Standard_deviation_and_∂standard_deviation∂parameter ∈ 4-element Vector{Symbol} And data, 2×4 Matrix{Float64}: (:Standard_deviation) (:σ) (:α) (:std_a) (:Pi) 0.0134547 -0.00218051 -0.0124107 0.561568 (:W_real) 0.174772 -0.0802741 0.193677 1.40047
for both inflation and real wages the volatility is higher and the standard deviation of the total factor productivity shock std_a has a much larger impact on the standard deviation of real wages compared to the linear solution.
At third order we get the following results:
julia> get_std(Gali_2015, parameter_derivatives = [:σ, :std_a, :α], variables = [:W_real,:Pi], algorithm = :pruned_third_order)Take symbolic derivatives up to third order: 1.902 seconds 2-dimensional KeyedArray(NamedDimsArray(...)) with keys: ↓ Variables ∈ 2-element Vector{Symbol} → Standard_deviation_and_∂standard_deviation∂parameter ∈ 4-element Vector{Symbol} And data, 2×4 Matrix{Float64}: (:Standard_deviation) (:σ) (:α) (:std_a) (:Pi) 0.0307732 0.00118823 0.0667707 3.55331 (:W_real) 0.260973 -0.148424 0.6022 9.27574
standard deviations of inflation is more than two times as high and for real wages it is also substantially higher. Furthermore, standard deviations of shocks matter even more for the volatility of the endogenous variables.
These results make it clear that capturing the nonlinear interactions by using nonlinear solutions has important implications for the model moments and by extension the model dynamics.
Method of moments for nonlinear solutions
Matching the theoretical moments of the nonlinear model solution to the data is no more complicated for the user than in the linear solution case (see above).
We need to define the target value and function and let an optimiser find the parameters minimising the distance to the target.
Keeping the targets:
- Mean of
W_real = 0.7 - Standard deviation of
Pi = 0.01
we need to define the target function and specify that we use a nonlinear solution algorithm (e.g. pruned third order):
julia> function distance_to_target(parameter_value_inputs) model_statistics = get_statistics(Gali_2015, parameter_value_inputs, algorithm = :pruned_third_order, parameters = [:α, :std_a], mean = [:W_real], standard_deviation = [:Pi]) targets = [0.7, 0.01] return sum(abs2, vcat(model_statistics...) - targets) enddistance_to_target (generic function with 1 method)
and then we can use the same code to optimise as in the linear solution case:
julia> sol = Optim.optimize(distance_to_target, [0,0], [1,1], [0.25, 0.01], Optim.Fminbox(Optim.LBFGS(linesearch = LineSearches.BackTracking(order = 3))))* Status: success * Candidate solution Final objective value: 2.703206e-05 * Found with Algorithm: Fminbox with L-BFGS * Convergence measures |x - x'| = 2.91e-06 ≰ 0.0e+00 |x - x'|/|x'| = 1.47e-05 ≰ 0.0e+00 |f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00 |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00 |g(x)| = 2.92e-09 ≤ 1.0e-08 * Work counters Seconds run: 15 (vs limit Inf) Iterations: 3 f(x) calls: 81 ∇f(x) calls: 38
the calculations take substantially longer and we don't get as close to our target as for the linear solution case. The parameter values minimising the distance are:
julia> sol.minimizer2-element Vector{Float64}: 0.19722918357256736 2.92352736050945e-9
lower than for the linear solution case and the theoretical moments given these parameter are:
julia> get_statistics(Gali_2015, sol.minimizer, algorithm = :pruned_third_order, parameters = [:α, :std_a], mean = [:W_real], standard_deviation = [:Pi])2-element Vector{AbstractArray{Float64}}: [0.6999560170816536] [0.015199050418199647]
The solution does not match the standard deviation of inflation very well.
Potentially the partial derivatives change a lot for small changes in parameters and even though the partial derivatives for standard deviation of inflation were large wrt std_a they might be small for values returned from the optimisation. We can check this with:
julia> get_std(Gali_2015, parameter_derivatives = [:σ, :std_a, :α], variables = [:W_real,:Pi], algorithm = :pruned_third_order, parameters = [:α, :std_a] .=> sol.minimizer)2-dimensional KeyedArray(NamedDimsArray(...)) with keys: ↓ Variables ∈ 2-element Vector{Symbol} → Standard_deviation_and_∂standard_deviation∂parameter ∈ 4-element Vector{Symbol} And data, 2×4 Matrix{Float64}: (:Standard_deviation) (:σ) (:α) (:std_a) (:Pi) 0.0151991 -0.00757065 -0.00698981 0.106398 (:W_real) 0.202943 -0.137097 0.313429 0.819094
and indeed it seems also the second derivative is large since the first derivative changed significantly.
Another parameter we can try is σ. It has a positive impact on the mean of real wages and a negative impact on standard deviation of inflation.
We need to redefine our target function and optimise it. Note that the previous call made a permanent change of parameters (as do all calls where parameters are explicitly set) and now std_a is set to 2.91e-9 and no longer 0.01.
julia> function distance_to_target(parameter_value_inputs) model_statistics = get_statistics(Gali_2015, parameter_value_inputs, algorithm = :pruned_third_order, parameters = [:α, :σ], mean = [:W_real], standard_deviation = [:Pi]) targets = [0.7, 0.01] return sum(abs2, vcat(model_statistics...) - targets) enddistance_to_target (generic function with 1 method)julia> sol = Optim.optimize(distance_to_target, [0,0], [1,3], [0.25, 1], Optim.Fminbox(Optim.LBFGS(linesearch = LineSearches.BackTracking(order = 3))))* Status: success * Candidate solution Final objective value: 1.567192e-12 * Found with Algorithm: Fminbox with L-BFGS * Convergence measures |x - x'| = 2.85e-08 ≰ 0.0e+00 |x - x'|/|x'| = 1.33e-08 ≰ 0.0e+00 |f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00 |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00 |g(x)| = 6.91e-09 ≤ 1.0e-08 * Work counters Seconds run: 13 (vs limit Inf) Iterations: 3 f(x) calls: 41 ∇f(x) calls: 41julia> sol.minimizer2-element Vector{Float64}: 0.20874375760768965 2.12374177099488
Given the new value for std_a and optimising over σ allows us to match the target exactly.