Comparing Models

This document will compare the solutions from NLP_model and MPSGE_model. They both describe the same economic model, but one is implemented in the traditional NLP framework and the other in MPSGE.

using NLP_to_MPSGE_Example
using JuMP, MPSGE, DataFrames

Initialization

Initialize the data

data = ModelData()

ModelData(25.0, 25.0, 100.0, 100.0, 75.0, 75.0, 1.0, 1.0, 1.0, 1.0)

Initialize the models

nlp = NLP_model(data);
mpsge = MPSGE_model(data);

Verify the models by solving the baseline scenario. For the NLP model, you are looking for Optimal Solution Found.

optimize!(nlp)

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.1.

Number of nonzeros in equality constraint Jacobian...:       48
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:       32

Total number of variables............................:       14
                     variables with only lower bounds:       12
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:       14
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.0000000e+02 7.46e-14 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.0000000e+02 1.42e-14 2.31e-15  -1.0 2.89e-14    -  9.91e-01 1.00e+00h  1

Number of Iterations....: 1

                                   (scaled)                 (unscaled)
Objective...............:  -1.0000000000000003e+02    1.0000000000000003e+02
Dual infeasibility......:   2.3093924602690858e-15    2.3093924602690858e-15
Constraint violation....:   1.4210854715202004e-14    1.4210854715202004e-14
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   1.4210854715202004e-14    1.4210854715202004e-14


Number of objective function evaluations             = 2
Number of objective gradient evaluations             = 2
Number of equality constraint evaluations            = 2
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 1
Total seconds in IPOPT                               = 3.998

EXIT: Optimal Solution Found.

For the MPSGE model, you are looking for Solver Status: LOCALL_SOLVED, or a Postsolved Residual that is very small (ideally, 0).

solve!(mpsge; cumulative_iteration_limit=0)

Reading options file /tmp/jl_AAG34b
 > cumulative_iteration_limit 0
Read of options file complete.

Path 5.0.03 (Fri Jun 26 09:54:13 2020)
Written by Todd Munson, Steven Dirkse, Youngdae Kim, and Michael Ferris
Preprocessed size   : 10

Major Iteration Log
major minor  func  grad  residual    step  type prox    inorm  (label)
    0     0     1     1 0.0000e+00           I 0.0e+00 0.0e+00 (z_p[X)

Major Iterations. . . . 0
Minor Iterations. . . . 0
Restarts. . . . . . . . 0
Crash Iterations. . . . 0
Gradient Steps. . . . . 0
Function Evaluations. . 1
Gradient Evaluations. . 1
Basis Time. . . . . . . 0.000002
Total Time. . . . . . . 0.240259
Residual. . . . . . . . 0.000000e+00
Postsolved residual: 0.0000e+00


Solver Status: LOCALLY_SOLVED
Model Status: FEASIBLE_POINT

At this point I'll set both models to silent mode to avoid cluttering the output

set_silent(nlp)
set_silent(mpsge)

Comparing Solutions

The NLP version of the model has more variables the MPSGE version, however, the solutions to both are the same.

Set parameter values for a counterfactual

params = ModelParameters(
    balance_of_payments = 10.0,
    elas_substitution = 0.5,
    elas_transformation = 0.5,
    tax_domestic = .1,
    tax_import = .2,
    price_world_export = 1.1,
    price_world_import = 1.2,
    subsidy_export = 0.05,
    )

set_parameter_values(nlp, params)
set_parameter_values(mpsge, params)

optimize!(nlp)
solve!(mpsge)

This is a method to view all the variable values in a DataFrame. Let's use this to compare the two models. Note that I've used some fancy operations, I've tried to space things out for clarity, you can extract each piece and run it separately to understand what is happening.

outerjoin(

    zip(JuMP.name.(all_variables(nlp)), value.(all_variables(nlp))) |>
        x -> DataFrame(var = first.(x), nlp = last.(x)),

    generate_report(mpsge) |>
        x -> transform(x, :var => ByRow(y -> JuMP.name(y)) => :var) |>
        x -> select(x, :var, :value => :mpsge),

    on = :var,
)

The variables present in the MPSGE model can be directly compared to the NLP model. However, neither X nor Q match. This is because in MPSGE, these are activity levels that produce quantities of goods, whereas in the NLP model, these are quantities of goods themselves. Let's focus on X as Q is similar.

Look at the MPSGE model and search for any productions that use X0. That would be the X block

production(mpsge[:X])

$Production: X
:t = omega
  O:PFX    Q:25 PWE    P:1.0 / PWE    A:Y    T:-TE
  O:PDD    Q:75.0

:s = 0
  I:PX    Q:100.0

Or, in the code:

@production(MP, X, [t=omega,s=0], begin @output(PFX, E0*PWE, t, reference_price = 1/PWE, taxes = [Tax(Y, -TE)]) # Negative tax for export subsidy @output(PDD, DS0, t) @input(PX, X0, s) end)

PX is the commodity associated with X0, so to get X we compute the compensated demand for X with respect to PX, and scale it by the activity level of X:

X = value(compensated_demand(mpsge[:X], mpsge[:PX])*mpsge[:X])
Q = value(compensated_demand(mpsge[:Q], mpsge[:PQ])*mpsge[:Q])

println("X = $X\nQ = $Q")

X = 100.0
Q = -105.87765111270654

The - on Q is signifying that the value comes from an output.

Next up, parameters. MPSGE does not report the parameter values. We can display these using:

value.(parameters(mpsge))

8-element Vector{Float64}:
  1.1
  0.1
  0.5
 10.0
  0.2
  0.05
  1.2
  0.5

The above displays only values. It would be a good exercise to put these in a DataFrame with the parameter names, and vcat it to the previous mpsge dataframe.

Next, the four quantities, DD, DS, E, and M. The idea is identical to X and Q. The exception is for E and M, which require scaling by the world prices.

DD = value(compensated_demand(mpsge[:Q], mpsge[:PDD])*mpsge[:Q])
DS = value(compensated_demand(mpsge[:X], mpsge[:PDD])*mpsge[:X])
E = value(compensated_demand(mpsge[:X], mpsge[:PFX])*mpsge[:X]/mpsge[:PWE])
M = value(compensated_demand(mpsge[:Q], mpsge[:PFX])*mpsge[:Q]/mpsge[:PWM])

println("DD = $DD\nDS = $DS\nE = $E\nM = $M")

DD = 77.5546334628556
DS = -77.55463346285558
E = -22.015179618694756
M = 28.513914650470234

The GR variable in the NLP model corresponds to the tax revenues in MPSGE. There should be a better function for this, but for now we can compute it as:

value(MPSGE.tax_revenue(mpsge[:X], mpsge[:Y]; virtual = true)) +
value(MPSGE.tax_revenue(mpsge[:Q], mpsge[:Y]; virtual = true))

-10.762120334674897

Finally, the three prices. I am first going to demonstrate a method to extract the price directly from the MPSGE model, then show how to compute them manually. This should be a function in MPSGE, I will work on it for a future release. However, notice that we can extract the expression, not just the value.

P = production(mpsge[:Q])
O = input(P)
N = O.children[2]
CF = cost_function(N)
PDT = value(CF)

println("PDT = $PDT\nCF = $CF")

PDT = 1.0483725073094727
CF = PDD*TD + PDD

Manually:

PDT = value(mpsge[:PDD]*(1 + mpsge[:TD]))
PMD = value(mpsge[:PFX]*(1 + mpsge[:TM])*mpsge[:PWM])
PED = value(mpsge[:PFX]*(1 + mpsge[:TE])*mpsge[:PWE])

println("PDT = $PDT\nPMD = $PMD\nPED = $PED")

PDT = 1.0483725073094727
PMD = 0.8617373614597765
PED = 0.6911851753375292

Recreating Results

In this part we will recreate tables 5.4 and 5.5 from the reference paper. These tables are on pages 42 (45 of document) and 45 (48 of document). These tables show the effects of changing the elasticity parameters and balance of payments on key variables in the model.

shocks = [
    ModelParameters(balance_of_payments = 0,    elas_substitution = 0.2, elas_transformation = 0.2,),
    ModelParameters(balance_of_payments = 10.0, elas_substitution = 0.2, elas_transformation = 0.2,),
    ModelParameters(balance_of_payments = 10.0, elas_substitution = 0.5, elas_transformation = 0.5,),
    ModelParameters(balance_of_payments = 10.0, elas_substitution = 2,   elas_transformation = 2,  ),
    ModelParameters(balance_of_payments = 10.0, elas_substitution = 5,   elas_transformation = 5,  ),
    ModelParameters(balance_of_payments = 10.0, elas_substitution = 15,  elas_transformation = .2, ),
    ModelParameters(balance_of_payments = 10.0, elas_substitution = .2,  elas_transformation = 15, ),

    ModelParameters(price_world_import = 1.1, elas_substitution = 0.2, elas_transformation = 0.2,),
    ModelParameters(price_world_import = 1.1, elas_substitution = 0.5, elas_transformation = 0.5,),
    ModelParameters(price_world_import = 1.1, elas_substitution = 2,   elas_transformation = 2,  ),
    ModelParameters(price_world_import = 1.1, elas_substitution = 5,   elas_transformation = 5,  ),
    ModelParameters(price_world_import = 1.1, elas_substitution = 15,  elas_transformation = .2, ),
    ModelParameters(price_world_import = 1.1, elas_substitution = .2,  elas_transformation = 15, ),
];

To run the shocks and collect the results, we do two steps:

1. Create empty DataFrames to hold the results
2. Loop through each shock, set the parameters, solve the models, generate the reports, and append them to the results DataFrames.

However, this requires the DataFrame be empty at the start of each loop. To ensure this, we can use a begin...end block to create a code block.

begin
    nlp_results = DataFrame();
    mpsge_results = DataFrame();
    for shock in shocks
        set_parameter_values(nlp, shock)
        set_parameter_values(mpsge, shock)

        optimize!(nlp)
        solve!(mpsge)

        nlp_report = report(nlp)
        mpsge_report = report(mpsge)

        push!(nlp_results, nlp_report[1, :])
        push!(mpsge_results, mpsge_report[1, :])
    end
end

The table results can be viewed by simply typing the DataFrame name.

nlp_results

MPSGE results

mpsge_results

Are these the same? Substract and round:

round.(nlp_results .- mpsge_results, digits=6)

What do the columns in the report mean? Check the doc string for report. You can do this either by entering the help mode in the REPL, just type ? in the REPL to switch to help mode and type report, or using the @doc report. Here is a link to the report documentation.

This page was generated using Literate.jl.