juliaacademy / decision-making-under-uncertainty Goto Github PK

# Parameters for the Neoclassical Growth Model (NGM)
α = 0.65; 
f(s;α=α) = (s)^α
a0(s)    = f(s)
min_s = eps(); max_s = 2.0; n_s = 100;
min_a = eps(); max_a = a0(max_s); n_a = 100; 

"states:"
states = range(min_s, max_s, length=n_s)
ceil_state(s) = states[searchsortedfirst(states, s)] # for discrete VFI

"actions:" 
valid_actions()  = range(min_a, max_a, length=n_a)        # possible actions any state.
valid_actions(s) = filter(<(a0(s)), valid_actions())      # valid actions @ state=s 

"Transition:" 
μ(s,a)      = f(s) - a      

"reward:" 
r(s,a) = log(a)

"discount:"
β = 0.95

using QuickPOMDPs: QuickMDP           #QuickMDP()
using POMDPModelTools: Deterministic
m = QuickMDP(
    states     = states,
    actions    = valid_actions,
    transition = (s, a) -> Deterministic( ceil_state( μ(s,a) ) ),
    reward     = (s, a) -> r(s,a),
    discount   = β
)

# DiscreteValueIteration: Both work fast enough
using DiscreteValueIteration
s1 = DiscreteValueIteration.SparseValueIterationSolver()
s2 = DiscreteValueIteration.ValueIterationSolver()
@time sol1 = solve(s1, m)
@time sol2 = solve(s2, m)
value(sol1, states[2]), action(sol1, states[2])
value(sol2, states[2]), action(sol2, states[2])

#rewrite MDP w/o ceil_state() in transition.
m = QuickMDP(
    states     = states,
    actions    = valid_actions,
    transition = (s, a) -> Deterministic( μ(s,a) ),
    reward     = (s, a) -> r(s,a),
    discount   = β
)

# LocalApproximationValueIteration works, but currently slow!
using GridInterpolations
using LocalFunctionApproximation
using LocalApproximationValueIteration
grid = GridInterpolations.RectangleGrid(states, valid_actions(), ) #[0.0, 1.0]) 
interp = LocalFunctionApproximation.LocalGIFunctionApproximator(grid)
s4 = LocalApproximationValueIterationSolver(interp)
@time sol4 = solve(s4, m)
#190.477798 seconds (2.37 G allocations: 82.610 GiB, 6.87% gc time, 0.21% compilation time)
value(sol4, states[2]), action(sol4, states[2])

using Plots
ω = α*β
A1 = α/(1.0-ω)
A0 = ((1-ω)*log(1-ω) + ω*log(ω))/((1-ω)*(1-β))

# Value
plot(legend=:bottomright, title="Value Functions");
plot!(states[2:end], i->A0 + A1 * log(i), lab="closed form") # way off 
plot!(states[2:end], i->value(sol1, i),   lab="sol1")
plot!(states[2:end], i->value(sol2, i),   lab="sol2")
plot!(states[2:end], i->value(sol4, i),   lab="sol4")
# Policy
plot(legend=:bottomright, title="Policy Functions");
plot!(states[2:end], i -> (1-ω)*(i^α), lab="closed form") # way off 
plot!(states[2:end], i -> action(sol1, i),   lab="sol1")
plot!(states[2:end], i -> action(sol2, i),   lab="sol2")
plot!(states[2:end], i -> action(sol4, i),   lab="sol4")
# Simulation
Tsim=150; s0=states[2]; 

simcf = []; push!(simcf, s0); 
for tt in 1:Tsim
    tt==1 ? s = s0 : nothing 
    a = (1-ω)*(s^α)
    sp = μ(s,a) 
    #sp = ceil_state(sp)
    push!(simcf, sp)
    s = sp
end 


sim1 = []; push!(sim1, s0); 
for tt in 1:Tsim
    tt==1 ? s = s0 : nothing 
    #a = valid_actions()[sol1.policy[searchsortedfirst(states, s)]]
    a = action(sol1, s)
    sp = μ(s,a) 
    sp = ceil_state(sp)
    push!(sim1, sp)
    s = sp
end 

sim4 = []; push!(sim4, s0); 
for tt in 1:Tsim
    tt==1 ? s = s0 : nothing 
    a = action(sol4, s)
    sp = μ(s,a) 
    push!(sim4, sp)
    s = sp
end 

plot(legend=:bottomright, title="Simulation");
plot!(simcf,   lab="closed form")
plot!(sim1,   lab="sol1")
plot!(sim4,   lab="sol4")