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 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, the package is loaded and then the @model macro is used to define the 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 are written 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 the expression is wrapped 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 the parameters of the model need to be added. 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.245 seconds Set up non-stochastic steady state problem: 0.586 seconds Find non-stochastic steady state: 1.258 seconds Take symbolic derivatives up to first order: 0.131 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 values are assigned (e.g. α = .25).
Note that one parameter definition per line is required.
Linear solution
Inspect model moments
Given the equations and parameters, everything is available for the package to generate the theoretical model moments. The mean of the linearised model can be retrieved as follows:
julia> get_mean(Gali_2015)[ Info: Most of the time is spent calculating derivatives wrt parameters. If they are not needed, add `derivatives = false` as an argument to the function call. 2-dimensional KeyedArray(NamedDimsArray(...)) with keys: ↓ Variables ∈ 23-element Vector{Symbol} → Mean_and_∂mean∂parameter ∈ 17-element Vector{Symbol} And data, 23×17 Matrix{Float64}: (:Mean) … (:std_a) (:std_z) (:std_nu) (:A) 1.0 0.0 0.0 0.0 (:C) 0.95058 0.0 0.0 0.0 (:MC) 0.888889 0.0 0.0 0.0 (:M_real) 0.915236 0.0 0.0 0.0 (:N) 0.934655 … 0.0 0.0 0.0 (:Pi) 1.0 0.0 0.0 0.0 ⋮ ⋱ ⋮ (:log_y) -0.0506831 0.0 0.0 0.0 (:nu) 0.0 0.0 0.0 0.0 (:pi_ann) -3.10862e-15 … 0.0 0.0 0.0 (:r_real_ann) 0.0402013 0.0 0.0 0.0 (:realinterest) 1.0101 0.0 0.0 0.0 (:x_aux_1) 3.452 0.0 0.0 0.0 (:x_aux_2) 3.8835 0.0 0.0 0.0
and the standard deviation like this:
julia> get_standard_deviation(Gali_2015)[ Info: Most of the time is spent calculating derivatives wrt parameters. If they are not needed, add `derivatives = false` as an argument to the function call. 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 1.91562e-30 3.18576e-32 (: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 1.38853e-29 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
Alternatively, std or get_std can be used to achieve 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_periods ∈ 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) 3.885e-18 -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 calibrating the model, examine how parameter changes affect model moments. MacroModelling.jl provides partial derivatives of model moments with respect to model parameters. This model is medium-sized, and derivatives are shown automatically. In this example, the sensitivity of the mean of all variables with respect to the production function parameter σ can be obtained like this:
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 -8.2121e-18 (:C) 0.95058 0.0060223 (:MC) 0.888889 -1.93054e-17 (:M_real) 0.915236 0.00579838 (:N) 0.934655 0.00789521 (:Pi) 1.0 3.16441e-19 ⋮ (:log_y) -0.0506831 0.00633539 (:nu) 0.0 0.0 (:pi_ann) -3.10862e-15 1.73472e-18 (:r_real_ann) 0.0402013 2.54586e-19 (:realinterest) 1.0101 6.42894e-20 (: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 -8.2121e-18 -8.89202e-16 1.89791e-16 (:C) 0.95058 0.0060223 3.79928e-15 -0.0941921 (:MC) 0.888889 -1.93054e-17 3.04393e-14 -1.28399e-15 (:M_real) 0.915236 0.00579838 3.48529 -0.09069 (:N) 0.934655 0.00789521 … 5.25019e-15 -0.207701 (:Pi) 1.0 3.16441e-19 4.86053e-16 6.11523e-34 ⋮ ⋱ ⋮ (:log_y) -0.0506831 0.00633539 3.88578e-15 -0.0990891 (:nu) 0.0 0.0 2.19824e-16 -9.53335e-35 (:pi_ann) -3.10862e-15 1.73472e-18 … 1.88738e-15 3.08149e-33 (:r_real_ann) 0.0402013 2.54586e-19 -4.0404 3.7821e-34 (:realinterest) 1.0101 6.42894e-20 -1.0203 9.55076e-35 (:x_aux_1) 3.452 0.174958 10.0544 -1.2868e-13 (:x_aux_2) 3.8835 0.196828 11.3112 1.02555e-16
The same can be done for standard deviation or variance, and all parameters:
julia> get_std(Gali_2015, parameter_derivatives = get_parameters(Gali_2015))[ Info: Most of the time is spent calculating derivatives wrt parameters. If they are not needed, add `derivatives = false` as an argument to the function call. 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 1.91562e-30 3.18576e-32 (: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 1.38853e-29 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))[ Info: Most of the time is spent calculating derivatives wrt parameters. If they are not needed, add `derivatives = false` as an argument to the function call. 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 -1.84565e-19 8.78945e-32 1.46173e-33 (: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 -4.06576e-21 8.01669e-32 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
This information can be used to calibrate certain values to targets. For example, assuming higher real wages (:W_real), and lower inflation volatility are desired. Since there are too many variables and parameters to be shown here, only a subset of them is printed:
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 3.16441e-19 6.11523e-34 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 it can be seen that lowering the production function parameter :α will increase real wages, but at the same time it will increase inflation volatility. This effect could be compensated by decreasing the standard deviation of the total factor productivity shock :std_a.
Method of moments
Instead of doing this by hand a target can also be set and an optimiser can find the corresponding parameter values. In order to do that targets need to be defined, and an optimisation problem needs to be set up.
The targets are:
- Mean of
W_real = 0.7 - Standard deviation of
Pi = 0.01
For the optimisation problem the L-BFGS algorithm implemented in Optim.jl is used. 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 get_statistics can be used 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])Dict{Symbol, AbstractArray{Float64}} with 2 entries: :mean => [0.678025] :standard_deviation => [0.0123588]
First the model object is passed on, followed by the parameter values and the parameter names the values correspond to. Then the desired outputs are defined: for the mean real wages are wanted and for the standard deviation inflation is wanted. Outputs for variance, covariance, or autocorrelation can also be obtained the same way as for the mean and standard deviation.
Next, a function measuring how close the model is to the target for given values of :α and :std_a can be defined:
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[:mean], model_statistics[:standard_deviation]) - targets) enddistance_to_target (generic function with 1 method)
Now the function can be tested with the current parameter values. In case the parameter values are not known they can also be looked 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
this allows testing the distance function:
julia> distance_to_target([0.25, 0.01])0.0004884527635317929
Next 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'| = 1.70e-07 ≰ 0.0e+00 |x - x'|/|x'| = 7.62e-07 ≰ 0.0e+00 |f(x) - f(x')| = 9.84e-16 ≰ 0.0e+00 |f(x) - f(x')|/|f(x')| = 1.93e-09 ≰ 0.0e+00 |g(x)| = 5.36e-09 ≤ 1.0e-08 * Work counters Seconds run: 0 (vs limit Inf) Iterations: 5 f(x) calls: 38 ∇f(x) calls: 33
The first argument to the optimisation call is the function defined previously, followed by lower and upper bounds, the starting values, and finally the algorithm. For the algorithm Fminbox has to be added because bounds are present (optional) and the specific line search method is set to speed up convergence (recommended but optional).
The output shows that the optimisation almost perfectly matches the target and the values of the parameters found by the optimiser are:
julia> sol.minimizer2-element Vector{Float64}: 0.22330255707150548 1.6375279113103097e-9
slightly lower for both parameters (in line with the previous insights from the sensitivities).
Combine the method of moments with estimation by adding the distance to the target as a penalty term to the posterior log-likelihood.
Nonlinear solutions
Up to this point the linearised solution of the model was used. 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 the mean for the pruned second order solution can be obtained as follows:
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.391 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.00566881 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. The effect of volatility can be seen 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 pruned third order perturbation is used by construction but the standard deviation does. Consider 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 the results are:
julia> get_std(Gali_2015, parameter_derivatives = [:σ, :std_a, :α], variables = [:W_real,:Pi], algorithm = :pruned_third_order)Take symbolic derivatives up to third order: 3.372 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.0307731 0.00118826 0.0667708 3.55329 (:W_real) 0.260973 -0.148424 0.602201 9.27575
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).
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
the target function needs to specify that a nonlinear solution algorithm is used (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[:mean], model_statistics[:standard_deviation]) - targets) enddistance_to_target (generic function with 1 method)
and then 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.703124e-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')| = 3.16e-09 ≰ 0.0e+00 |f(x) - f(x')|/|f(x')| = 1.17e-04 ≰ 0.0e+00 |g(x)| = 2.92e-09 ≤ 1.0e-08 * Work counters Seconds run: 17 (vs limit Inf) Iterations: 3 f(x) calls: 51 ∇f(x) calls: 37
the calculations take substantially longer and the solution does not get as close to the target as for the linear solution case. The parameter values minimising the distance are:
julia> sol.minimizer2-element Vector{Float64}: 0.19722805576968225 2.924044828248831e-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])Dict{Symbol, AbstractArray{Float64}} with 2 entries: :mean => [0.699956] :standard_deviation => [0.015199]
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. This can be checked 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.015199 -0.00757059 -0.00699053 0.106387 (:W_real) 0.202941 -0.137095 0.31342 0.818991
and indeed it seems also the second derivative is large since the first derivative changed significantly.
Another parameter to try is σ. It has a positive impact on the mean of real wages and a negative impact on standard deviation of inflation.
The target function needs to be redefined and optimised. 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[:mean], model_statistics[:standard_deviation]) - 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.566886e-12 * Found with Algorithm: Fminbox with L-BFGS * Convergence measures |x - x'| = 2.55e-08 ≰ 0.0e+00 |x - x'|/|x'| = 1.20e-08 ≰ 0.0e+00 |f(x) - f(x')| = 1.25e-11 ≰ 0.0e+00 |f(x) - f(x')|/|f(x')| = 8.89e-01 ≰ 0.0e+00 |g(x)| = 6.91e-09 ≤ 1.0e-08 * Work counters Seconds run: 17 (vs limit Inf) Iterations: 3 f(x) calls: 41 ∇f(x) calls: 41julia> sol.minimizer2-element Vector{Float64}: 0.20874123864588165 2.1236961104956027
Given the new value for std_a and optimising over σ allows matching the target exactly.