Orbits/Code/2Firm_vfi.jl

#=
The purpose of this is to set up a 2-firm, finite, dynamic orbit use model.

# TODO LIST
 - [ ] Set up loss functions
 - [ ] Set up utility functions
 - [ ] Set up State×Choice → Utility matrix
 - [ ] Set up VFI
 - [ ] Run it
 - [ ] ?
 - [ ] Profit!!!!
=#

#=
Model Description

V({sʲₜ}, Dₜ) = max{xⁱ} U(sⁱₜ) - F(xⁱ) + β V({sʲₜ₊₁}, Dₜ₊₁)

Utility Function: Similar to Rao Rondina
U(s) = ηs - ρ

Laws of motion
sʲₜ₊₁ = (1-l({sʲₜ},Dₜ))sʲₜ + xʲₜ
Dₜ₊₁ = (1-δ)Dₜ + g(Dₜ) + γ⋅ ∑ⱼ sʲₜ⋅l({sʲₜ},Dₜ) + Γ⋅∑ xʲₜ

Loss Function: Similar to Rao Rondina
lⁱ({sʲₜ},Dₜ) = 1-exp(  -α₁⋅∑{i≠j}sʲₜ - α₂⋅sʲₜ - α₃⋅Dₜ  )

Autofragmentation Function: See Rao Rondina for a description (pg48)
g(Dₜ) = ζ⋅Dₜ⋅(1-exp(-α₄⋅Dₜ))

Launch Costs
F(xⁱ) = F⋅xⁱ

=#

#Setup Model Constants
const α1 = 1e-2
const α2 = 1e-2
const α3 = 1e-2
const α4 = 1e-2
const ζ = 3
const δ = 2e-2
const F = 5
const γ = 10
const Γ = 2

#=
Potential Structs

Maybe a view onto the state space?
With methods which show own and other sizes as vectors?
=#

#Loss Function
function loss(constellation_sizes::Vector{<:Int},debris::Float64,owner::Int)
    1-exp(-α1*sum(constellation_sizes) + (α1-α2)*constellation_sizes[owner] - α3*debris)
end

#=
I think that I might need to consider the launch quantities  of the other operators
as a state variable and use a nash equilibrium style idea of equilibrium.
=#