Simple Firm investment problem

[azev77]

(This is an edited/corrected version of my original post!)

I would like to contribute this example (once I confirm it works).

Consider a firm endowed with initial capital K0, that chooses investment to maximize the present value of its future dividends.
The firm terminates production at time T and sells its remaining undepreciated capital (or scraps it if its cheaper) .

The appropriate terminal condition for i[T] depends on the technology parameter z:

In principle, the terminal conditions can be combined into one:
@constraint(m, i == max(-1, -(1-δ)*k), DomainRestrictions(t => T))
This currently gives an error.

Consider the closed form for the case z≥z̄ or i[T]=-1:

Now w/ InfiniteOpt.jl

using InfiniteOpt, Ipopt, Plots
δ=0.10;  # capital depreciation
ρ=0.05;  # discount rate
T = 90.0 # life horizon
k0 = 1.0 # endowment
# Terminal Condition: Burn/Liquidate Capital. 
znum = (ρ+δ)*(ρ+2*δ)*exp((ρ+δ)*T)*(exp(δ*T)-(1-δ)*(1+δ*k0))
zden = (1-δ)*(δ + (ρ+δ)*exp((ρ+2*δ)*T) - (ρ+2*δ)*exp((ρ+δ)*T) )
zqot = znum/zden
#
z=zqot +0.05    
z > ρ + δ # Need: z > ρ + δ
I_SS = (z-ρ-δ)/(ρ+δ); 
K_SS = I_SS/δ;

# Closed Form Solution if `z ≥ zqot` or `i[T]=-1`:
ii(t) = I_SS - (1+I_SS)*exp(-(ρ+δ)*(T-t))
kk(t) = K_SS - 
    ((1+I_SS)/(ρ+2*δ))*exp(-(ρ+δ)*(T-t)) + 
    (((1+I_SS)/(ρ+2*δ))*exp(-(ρ+δ)*T) + k0 - K_SS)*exp(-(δ)*t);
#

u(k,i;z=z) = z*k -i -(1/2)*(i)^2 # utility function
discount(t; ρ=ρ) = exp(-ρ*t)     # discount function
BC(k, i; δ=δ) = i - δ*k          # LOM 

opt = Ipopt.Optimizer            # desired solver
ns = 10_000;                     # number of gridpoints

m = InfiniteModel(opt)
@infinite_parameter(m, t ∈ [0, T], num_supports = ns)
@variable(m, k, Infinite(t)) ## state variables
@variable(m, i, Infinite(t)) ## control variables
@objective(m, Max, integral(u(k,i), t, weight_func = discount))
@constraint(m, c1, deriv(k, t) == BC(k, i; δ=δ))
@constraint(m, c2, i >= -(1-δ)*k)                   # Can't sell more than all your capital.
@constraint(m, k == k0, DomainRestrictions(t => 0))
@constraint(m, i >=-1, DomainRestrictions(t => T))
@constraint(m, i >=-(1-δ)*k, DomainRestrictions(t => T)) # redundant w/ c2 

optimize!(m)
termination_status(m)

i_opt = value(i)
k_opt = value(k)
ts = supports(t)
opt_obj = objective_value(m) # V(k0, 0)

ix = 2:(length(ts)-1) # index for plotting
#
plot(legend=:topright);
plot!(ts[ix], i_opt[ix], color = 1, lab = "i: InfiniteOpt")
plot!(ts[ix], ii,        color = 1, linestyle=:dash, lab = "i: closed form")
plot!(ts[ix], k_opt[ix], color = 4, lab = "k: InfiniteOpt");
plot!(ts[ix], kk,        color = 4, linestyle=:dash, lab = "k: closed form")
plot!([0.0], seriestype = :hline, lab="", color="red")
plot!(ts[ix],   -(1-δ)*k_opt[ix], lab = "i>-(1-δ)*k", color="black", l=:dash)
plot!([-1],  seriestype = :hline, lab="-1",           color="black", l=:dash)

The solution in this case is very very accurate!

Now solve the case z<z̄ or i[T]=-(1-δ)*k[T]>-1:

using InfiniteOpt, Ipopt, Plots
δ=0.10;  # capital depreciation
ρ=0.05;  # discount rate
T = 90.0 # life horizon
k0 = 1.0 # endowment
# Terminal Condition: Burn/Liquidate Capital. 
znum = (ρ+δ)*(ρ+2*δ)*exp((ρ+δ)*T)*(exp(δ*T)-(1-δ)*(1+δ*k0))
zden = (1-δ)*(δ + (ρ+δ)*exp((ρ+2*δ)*T) - (ρ+2*δ)*exp((ρ+δ)*T) )
zqot = znum/zden
#
z=zqot - 0.10    
z > ρ + δ # Need: z > ρ + δ
I_SS = (z-ρ-δ)/(ρ+δ); 
K_SS = I_SS/δ;

# Closed Form Solution if `z ≥ zqot` or `i[T]=-1`:
ii(t) = I_SS - (1+I_SS)*exp(-(ρ+δ)*(T-t))
kk(t) = K_SS - 
    ((1+I_SS)/(ρ+2*δ))*exp(-(ρ+δ)*(T-t)) + 
    (((1+I_SS)/(ρ+2*δ))*exp(-(ρ+δ)*T) + k0 - K_SS)*exp(-(δ)*t);
#

u(k,i;z=z) = z*k -i -(1/2)*(i)^2 # utility function
discount(t; ρ=ρ) = exp(-ρ*t)     # discount function
BC(k, i; δ=δ) = i - δ*k          # LOM 

opt = Ipopt.Optimizer;           # desired solver
ns = 10_000;                     # number of gridpoints

m = InfiniteModel(opt)
@infinite_parameter(m, t ∈ [0, T], num_supports = ns)
@variable(m, k, Infinite(t)) ## state variables
@variable(m, i, Infinite(t)) ## control variables
@objective(m, Max, integral(u(k,i), t, weight_func = discount))
@constraint(m, c1, deriv(k, t) == BC(k, i; δ=δ))
@constraint(m, c2, i >= -(1-δ)*k)                   # Can't sell more than all your capital.
@constraint(m, k == k0, DomainRestrictions(t => 0))
@constraint(m, i >=-1, DomainRestrictions(t => T))
@constraint(m, i >=-(1-δ)*k, DomainRestrictions(t => T)) # redundant w/ c2 

optimize!(m)
termination_status(m)

i_opt = value(i)
k_opt = value(k)
ts = supports(t)
opt_obj = objective_value(m) # V(k0, 0)

ix = 2:(length(ts)-1) # index for plotting
#
plot(legend=:topright);
plot!(ts[ix], i_opt[ix], color = 1, lab = "i: InfiniteOpt")
# plot!(ts[ix], ii,        color = 1, linestyle=:dash, lab = "i: closed form")
plot!(ts[ix], k_opt[ix], color = 4, lab = "k: InfiniteOpt")
# plot!(ts[ix], kk,        color = 4, linestyle=:dash, lab = "k: closed form")
plot!([0.0], seriestype = :hline, lab="", color="red")
plot!(ts[ix],   -(1-δ)*k_opt[ix], lab = "i>-(1-δ)*k", color="black", l=:dash)
plot!([-1],  seriestype = :hline, lab="-1",           color="black", l=:dash)

Now the solution (for investment) has a weird unexpected kink:

Now plot it w/ the "i>-(1-δ)*k" constraint:

Compare this to the closed form solution in Mathematica (no kink & the expected constraint binds at "T"):

Do you know what I might be doing wrong?
My guess is I need to tune one of the optimization options...

[azev77]

Also I attempted to solve an infinite horizon version using methods you proposed here.

# CASE: Infinite Horizon  
using InfiniteOpt, Ipopt, Plots
δ=0.10;  # capital depreciation
ρ=0.05;  # discount rate
#T = 90.0 # life horizon
k0 = 1.0 # endowment
#
z= ρ + δ +0.10    
z > ρ + δ # Need: z > ρ + δ
I_SS = (z-ρ-δ)/(ρ+δ); 
K_SS = I_SS/δ;

# Closed Form Solution:
ii(t) = I_SS 
kk(t) = K_SS + exp(-(δ)*t)*(k0 - K_SS);
#

u(k,i;z=z) = z*k -i -(1/2)*(i)^2 # utility function
discount(t; ρ=ρ) = exp(-ρ*t)     # discount function
int_discount(t; ρ=ρ) = exp((1-ρ)*t) 
# additional weight function for the integral (will be multiplied by e^-t by Gauss Laguerre)
BC(k, i; δ=δ) = i - δ*k          # LOM 

opt = Ipopt.Optimizer            # desired solver
ns = 10_000;                     # number of gridpoints

m = InfiniteModel(opt)
@infinite_parameter(m, t in [0, Inf])
#@infinite_parameter(m, t ∈ [0, T], num_supports = ns)
@variable(m, k, Infinite(t)) ## state variables
@variable(m, i, Infinite(t)) ## control variables
@parameter_function(m, dis == discount(t))
@objective(m, Max, integral(u(k,i), t, eval_method = GaussLaguerre(), num_nodes = 100, weight_func = int_discount))
@constraint(m, k == k0, DomainRestrictions(t => 0))

@constraint(m, dis * (1+i) * k == 0, DomainRestrictions(t => supports(t)[end])) 
# not sure what u'(c_t) is...

@constraint(m, c1, deriv(k, t) == BC(k, i; δ=δ))
#@constraint(m, c2, i >= -(1-δ)*k)
optimize!(m)
termination_status(m)
opt_obj = objective_value(m) #
i_opt, k_opt = value(i), value(k)
ts = supports(t)

ix = 2:(length(ts)-1) # index for plotting
#
plot(legend=:bottomleft);
plot!(ts[ix], i_opt[ix], color = 1, lab = "i: InfiniteOpt")
plot!(ts[ix], ii,        color = 1, linestyle=:dash, lab = "i: closed form")
plot!(ts[ix], k_opt[ix], color = 4, lab = "k: InfiniteOpt");
plot!(ts[ix], kk,        color = 4, linestyle=:dash, lab = "k: closed form")
plot!([0.0], seriestype = :hline, lab="", color="red")
plot!(ts[ix],   -(1-δ)*k_opt[ix], lab = "i>-(1-δ)*k", color="black", l=:dash)
plot!([-1],  seriestype = :hline, lab="-1",           color="black", l=:dash)

But the solution isn't quite right:

[pulsipher]

Hi there!

First, good catch on the max bug. That seems to be a bug with JuMP, I have reached out to them to resolve it. So once that is fixed, your terminal constraint should work as expected.

Second, the "kink" behavior occurs because i(t) hits the constraint as you show. This all occurs because k(t) decays a little further relative to the simulation. However, I tried a variety of higher fidelity derivative approximations and get the same simulated response. Hence, the optimization problem is behaving as expected (in terms of respecting its constraints). Moreover, the problem is convex so this cannot be attributed to a local solution of Ipopt. So the simultaneous approach we use to solve the problem is behaving as I would expect, so I am not sure why there is a difference. There might exist a more problem specific reason.

Third, the behavior of the infinite horizon case isn't surprising since simultaneous approximation methods (which is currently what I have implemented) tend toward higher accuracy in the early times infinite horizon problems, but exhibit drift afterward because of the small weights you get in the integral (since the additional discount factor decreases the accuracy of the quadrature). More sophisticated solution methods are needed for infinite horizon problems, see #190.

Fourth, instead of using DomainRestrictions at particular points, you can use our restricted variable syntax.

@constraint(m, i >=-(1-δ)*k, DomainRestrictions(t => T))
@constraint(m, i(T) >=-(1-δ)*k(T)) # this is the same but will be more computationally efficient

[pulsipher]

Following up, #195 fixes the max bug for now. It must be registered under a different name until JuMP.jl releases its next version. Then it can just be registered as normal.

[azev77]

Hmm, I might be making a math mistake moving from discrete time to continuous time in the kink example.

In discrete time:
Recall the law of motion (LOM) for capital is: k[t+1] = i[t] +(1-δ)*k[t]
If we hit the liquidation constraint at time "j": i[j] =-(1-δ)*k[j]
If we plug that into the LOM then capital automatically goes to zero: k[j+1] =0

In continuous time:
Recall the law of motion (LOM) for capital is: k'[t] = i[t] -δ*k[t]
Briefly before time is over, we have i[j] =-(1-δ)*k[j]
If we plug that into the LOM: k'[j] = -k[j]
Then shouldn't capital immediately go to zero instantaneously after the i[t] >=-(1-δ)*k[t] constraint binds???

Perhaps the terminal condition in continuous time should instead be: k[T] =0
Note: in discrete time the terminal condition is k[T+1]=0
But when I solve w/ infiniteOpt.jl, the i[t] >=-(1-δ)*k[t] constraint appears to be violated.
hmmmm...

[azev77]

I think I figured out the problem in the case where all capital is liquidated.

in discrete time the terminal condition is: k_T+1 =0

continuous time analog: k_T=0
My closed form solution in this case assumed the i[t] >=-(1-δ)*k[t] constraint is slack. InfiniteOpt.jl knows it should bind briefly before the end of time..
Hence, the correct constrained solution in this case should still have a kink.

either way, economists typically only solve infinite horizon firm investment models…

[pulsipher]

Cool, that explains it then. The infinite horizon case is certainly of interest, but unfortunately I do not currently have the time needed to personally address #190. Contributions are more then welcome.

[azev77]

Here is the final answer to the finite-horizon firm investment problem.

The constrained optimization problem:

The technology threshold at which both terminal conditions bind:

Case 1: z = 0.30 > z̄ (burn capital)
Closed form solution:

InfiniteOpt.jl

using InfiniteOpt, Ipopt, Plots
δ=0.10;  # capital depreciation
ρ=0.05;  # discount rate
T = 90.0 # life horizon
k0 = 1.0 # endowment
z=0.3    
z > ρ + δ # Need: z > ρ + δ
I_SS = (z-ρ-δ)/(ρ+δ); 
K_SS = I_SS/δ;

# Closed Form Solution if `z ≥ zqot` or `i[T]=-1`:
ii(t) = I_SS - (1+I_SS)*exp(-(ρ+δ)*(T-t))
kk(t) = K_SS - 
    ((1+I_SS)/(ρ+2*δ))*exp(-(ρ+δ)*(T-t)) + 
    (((1+I_SS)/(ρ+2*δ))*exp(-(ρ+δ)*T) + k0 - K_SS)*exp(-(δ)*t);
#

u(k,i;z=z) = z*k -i -(1/2)*(i)^2 # utility function
discount(t; ρ=ρ) = exp(-ρ*t)     # discount function
BC(k, i; δ=δ) = i - δ*k          # LOM 

opt = Ipopt.Optimizer            # desired solver
ns = 10_000;                     # number of gridpoints

m = InfiniteModel(opt)
@infinite_parameter(m, t ∈ [0, T], num_supports = ns)
@variable(m, k, Infinite(t)) ## state variables
@variable(m, i, Infinite(t)) ## control variables
@objective(m, Max, integral(u(k,i), t, weight_func = discount))
@constraint(m, c1, deriv(k, t) == BC(k, i; δ=δ))
@constraint(m, c2, k >= 0)                   # Can't sell more than all your capital.
@constraint(m, k == k0, DomainRestrictions(t => 0))
@constraint(m, i(T) >=-1)
@constraint(m, k(T) >=0)


optimize!(m)
termination_status(m)

i_opt = value(i)
k_opt = value(k)
ts = supports(t)
opt_obj = objective_value(m) # V(k0, 0)

ix = 2:(length(ts)-1) # index for plotting
#
plot(legend=:topleft);
plot!(ts[ix], i_opt[ix], color = 1, lab = "i: InfiniteOpt")
plot!(ts[ix], ii,        color = 1, linestyle=:dash, lab = "i: closed form")
plot!([-1],  seriestype = :hline, lab="i>=-1, binds",      color=1, l=:dash)
plot!(ts[ix], k_opt[ix], color = 4, lab = "k: InfiniteOpt");
plot!(ts[ix], kk,        color = 4, linestyle=:dash, lab = "k: closed form")
plot!([0.0], seriestype = :hline, lab="k>=0, slack", color=4)

Case 2: z = 0.20 < z̄ (sell all capital at the end)
Closed form solution:

InfiniteOpt.jl

using InfiniteOpt, Ipopt, Plots
δ=0.10;  # capital depreciation
ρ=0.05;  # discount rate
T = 90.0 # life horizon
k0 = 1.0 # endowment
z=0.2    
z > ρ + δ # Need: z > ρ + δ
I_SS = (z-ρ-δ)/(ρ+δ); 
K_SS = I_SS/δ;

u(k,i;z=z) = z*k -i -(1/2)*(i)^2 # utility function
discount(t; ρ=ρ) = exp(-ρ*t)     # discount function
BC(k, i; δ=δ) = i - δ*k          # LOM 

opt = Ipopt.Optimizer            # desired solver
ns = 10_000;                     # number of gridpoints

m = InfiniteModel(opt)
@infinite_parameter(m, t ∈ [0, T], num_supports = ns)
@variable(m, k, Infinite(t)) ## state variables
@variable(m, i, Infinite(t)) ## control variables
@objective(m, Max, integral(u(k,i), t, weight_func = discount))
@constraint(m, c1, deriv(k, t) == BC(k, i; δ=δ))
@constraint(m, c2, k >= 0)                   # Can't sell more than all your capital.
@constraint(m, k == k0, DomainRestrictions(t => 0))
@constraint(m, i(T) >=-1)
@constraint(m, k(T) >=0)


optimize!(m)
termination_status(m)

i_opt = value(i)
k_opt = value(k)
ts = supports(t)
opt_obj = objective_value(m) # V(k0, 0)

ix = 2:(length(ts)-1) # index for plotting
#
plot(legend=:topleft);
plot!(ts[ix], i_opt[ix], color = 1, lab = "i: InfiniteOpt")
plot!([-1],  seriestype = :hline, lab="i>=-1, slack",      color=1, l=:dash)
plot!(ts[ix], k_opt[ix], color = 4, lab = "k: InfiniteOpt");
plot!([0.0], seriestype = :hline, lab="k>=0, binds", color=4)

PS: I don't know how to solve the infinite horizon version in your package.
I use many other methods in this discourse post.

At some point this can make for another nice example for the docs.

[azev77]

Here are some codes that solve the infinite horizon problem in continuous time
1: my HJB solver
2: using @matthieugomez EconPDEs.jl here
3: the ODEs can also be solved w/ your favorite DiffEq solver
I'm not sure if & to what extent this is helpful to InfiniteOpt.jl

[pulsipher]

Thanks for the info. The challenge here will be developing a solution method that can handle general infinite-dimensional optimization problems (not just time ones) with general differential algebraic equations and path/point constraints. Moreover, traditional ODE solvers are not typically amendable to optimization problems (due to the presence of decision variables), though they can be used with multiple shooting methods (which don't scale well unfortunately). This is why we currently implement simultaneous approaches to handle general DAEs. This works well for domains with finite bounds (e.g., finite time horizons), but extending it to domains with infinite bounds is nontrivial in general.