Comparative Statics w/ InfiniteOpt.jl

[azev77]

The Consumption Savings example is solved for one vector of parameters.
I want to do comparative statics: study how the solution changes as the parameter eg ρ changes.

My instinct is to define a function of the parameters, that solves the problem at a given set of parameters:

using InfiniteOpt, Ipopt

ρ = 0.025  # discount rate
k = 100.0  # utility bliss point
T = 10.0   # life horizon
r = 0.05   # interest rate
B0 = 100.0 # endowment
opt = Ipopt.Optimizer   # desired solver
ns = 1_000;             # number of gridpoints

function f(;ρ=ρ,k=k,T=T,r=r,B0=B0,opt=opt,ns=ns)
    u(c; k=k) = -(c - k)^2       # utility function
    discount(t; ρ=ρ) = exp(-ρ*t) # discount function
    BC(B, c; r=r) = r*B - c      # budget constraint
    #
    m = InfiniteModel(opt)
    @infinite_parameter(m, t in [0, T], num_supports = ns) 
    @variable(m, B, Infinite(t)) ## state variables
    @variable(m, c, Infinite(t)) ## control variables
    @objective(m, Max, integral(u(c), t, weight_func = discount))
    @constraint(m, B == B0, DomainRestrictions(t => 0))
    @constraint(m, B == 0, DomainRestrictions(t => T))
    @constraint(m, c1, deriv(B, t) == BC(B, c; r=r))
    #
    optimize!(m)
    termination_status(m)
    #
    c_opt = value(c)
    B_opt = value(B)
    #
    return [c_opt B_opt]    
end
# Define a grid of the parameter of interest ρ and solve over the grid:
grid_ρ = [.025, .05, .075]
res = []
for ρ in grid_ρ
    push!(res, f(ρ=ρ))
end 

This doesn't feel like a smart way to do comparative statics (don't know what they call it in your field).
Every time we solve the problem, we don't use the previous solution as a guess etc.
Is there a smarter/more powerful way to do this using your package?

I'm only asking in case there is a known, efficient, better approach to doing this.
Else don't waste time on this...

[pulsipher]

If I understand correctly you are interested in exploring the effect of ρ on the solution. This simply amounts to deriving a Pareto frontier, since each problem is fully decoupled from one another.

The code you have will work just fine, but I would probably just define the model once and update the objective function on the fly to better the performance:

using InfiniteOpt, Ipopt

# Define the parameters
k = 100.0  # utility bliss point
T = 10.0   # life horizon
r = 0.05   # interest rate
B0 = 100.0 # endowment
ns = 1_000;             # number of gridpoints

# Define expression functions
u(c; k=k) = -(c - k)^2       # utility function
BC(B, c; r=r) = r*B - c      # budget constraint

# Setup the model (without the objective)
m = InfiniteModel(Ipopt.Optimizer)
@infinite_parameter(m, t in [0, T], num_supports = ns) 
@variable(m, B, Infinite(t)) ## state variables
@variable(m, c, Infinite(t)) ## control variables
@constraint(m, B == B0, DomainRestrictions(t => 0))
@constraint(m, B == 0, DomainRestrictions(t => T))
@constraint(m, c1, deriv(B, t) == BC(B, c; r=r))

# Define a grid of the parameter of interest ρ and solve over the grid:
grid_ρ = [.025, .05, .075]
c_data = Dict()
B_data = Dict()
for ρ in grid_ρ
    @objective(m, Max, integral(u(c), t, weight_func = t -> exp(-ρ*t)))
    optimize!(m)
    c_data[ρ]= value(c)
    B_data[ρ]= value(B)
end 

Note that I haven't tested the above so it may have a typo somewhere.

[azev77]

Thanks!
Suppose I have a fine grid: grid_ρ = 0.025:0.001:0.075
From theory, I expect the solution at ρ=0.025 to be "close" to the solution at the next grid point ρ=0.026.

As I solve the model over this grid, is there a way to use the solution at grid_ρ[i-1] as a guess for the problem at grid_ρ[i-1]?

[pulsipher]

You can use the set_start_value_function method in the for loop to update the start (guess) values of the variables (see https://pulsipher.github.io/InfiniteOpt.jl/v0.4.0/guide/variable/#Start-Values). This can either be a constant value or a function (i.e., one that encodes the optimal trajectory). We don't however, currently provide a way to convert to the discretized solutions to functions (though that would make for a nice contribution and would relate to #82 and #44). However, you can manually accomplish this mapping using Interpolations.jl.

For example, we could use linear interpolation:

using InfiniteOpt, Ipopt, Interpolations

# Define the parameters
k = 100.0  # utility bliss point
T = 10.0   # life horizon
r = 0.05   # interest rate
B0 = 100.0 # endowment
ns = 1_000;             # number of gridpoints

# Define expression functions
u(c; k=k) = -(c - k)^2       # utility function
BC(B, c; r=r) = r*B - c      # budget constraint

# Setup the model (without the objective)
m = InfiniteModel(Ipopt.Optimizer)
@infinite_parameter(m, t in [0, T], num_supports = ns) 
@variable(m, B, Infinite(t)) ## state variables
@variable(m, c, Infinite(t)) ## control variables
@constraint(m, B == B0, DomainRestrictions(t => 0))
@constraint(m, B == 0, DomainRestrictions(t => T))
@constraint(m, c1, deriv(B, t) == BC(B, c; r=r))

# Define a grid of the parameter of interest ρ and solve over the grid:
grid_ρ = [.025, .05, .075]
ts = supports(t)
c_data = Dict()
B_data = Dict()
for ρ in grid_ρ
    # solve and save the results
    @objective(m, Max, integral(u(c), t, weight_func = t -> exp(-ρ*t)))
    optimize!(m)
    c_data[ρ]= value(c)
    B_data[ρ]= value(B)

    # setup up the guess values for the next iteration
    set_start_value_function(c, LinearInterpolation(ts, c_data[ρ]))
    set_start_value_function(B, LinearInterpolation(ts, B_data[ρ]))
end 

Again, I haven't tested this.

[azev77]

This line gives an error
set_start_value_function(c, LinearInterpolation(ts, c_data[ρ]))

julia>         set_start_value_function(c, LinearInterpolation(ts, c_data[ρ]))
ERROR: ArgumentError: `set_start_value_function` not defined for variable reference type `InfiniteVariableRef`.
Stacktrace:
 [1] set_start_value_function(vref::InfiniteVariableRef, start::Interpolations.Extrapolation{Float64, 1, Interpolations.GriddedInterpolation{Float64, 1, Float64, Gridded{Linear}, Tuple{Vector{Float64}}}, Gridded{Linear}, Throw{Nothing}})
   @ InfiniteOpt C:\Users\azevelev\.julia\packages\InfiniteOpt\8U5IK\src\general_variables.jl:718
 [2] set_start_value_function(vref::GeneralVariableRef, start::Interpolations.Extrapolation{Float64, 1, Interpolations.GriddedInterpolation{Float64, 1, Float64, Gridded{Linear}, Tuple{Vector{Float64}}}, Gridded{Linear}, Throw{Nothing}})
   @ InfiniteOpt C:\Users\azevelev\.julia\packages\InfiniteOpt\8U5IK\src\general_variables.jl:730
 [3] top-level scope
   @ REPL[73]:1

[pulsipher]

Oh I see, the problem is that LinearInterpolation doesn't give a function and that causes a method error. I will have to resolve that erroneous error you are getting there though.

You should be able to remedy the problem via a barrier function:

c_interp = LinearInterpolation(ts, c_data[ρ])
set_start_value_function(c, t -> c_interp(t))