Akira Matsushita
Department of Economics, University of Tokyo
From Miranda and Fackler, Applied Computational Economics and Finance, 2002, Section 8.4.8 and 9.7.8
Then the profit maximization problem for producer is formulated as:
\begin{align*} \underset{ \{x_t \geq 0 \}_{t=1}^T}{\max} \delta^T p s_{T+1} - \sum_{t=1}^{T} \delta^{t-1} \kappa x_t \hspace{1em} \text{s.t.} \hspace{1em} \begin{cases} s_{t+1} = g(s_t, x_t) \\ s_0 = \bar{s} \end{cases} \end{align*}where $\delta \geq 0$ is the discount factor
The corresponding Bellman equation is:
\begin{align*} \begin{cases} V_t(s_t) = \underset{x_t \geq 0}{\max} \{ -\kappa x_t + \delta V_{t+1}(g(s_t, x_t)) \} \hspace{1em} t=1, \ldots, T \\ V_{T+1}(s_{T+1}) = ps_{T+1} \end{cases} \end{align*}Let's derive the Euler Equilibrium Conditions from the Bellman equation
Assume $g: \mathbb{R}_+ \times \mathbb{R}_+ \rightarrow \mathbb{R}_+$ is differentiable on the entire domain (so $g$ is partially differentiable by both $s$ and $x$)
Define $g_x = \frac{\partial g}{\partial x}$ and $g_s = \frac{\partial g}{\partial s}$
Assume $g_x(s, 0)$ is large enough so that the nonnegative constraint of $x$ will not be binded at an optimum (ensures an interior point solution)
Let $X_t^*$ be the optimal solution correspondence, i.e.,
\begin{align*} X_t^*(s) = \left\{ x_t^* \in [0, \infty) \mid V_t(s) = -\kappa x_t^* + \delta V_{t+1}(g(s_t, x_t^*)) \right\} \end{align*}and $x_t^*(s) \in X_t^*(s)$ be a selection of it. Assume $x_t^*(s)$ is differentiable for all $t$
Note that by using $x_t^*(s)$ the optimal value function at time 1 is derived as
\begin{align*} V_1(\bar{s}) &= -\kappa x_1^*(\bar{s}) + \delta V_2(g(\bar{s}, x_1^*(\bar{s}))) \\ &= -\kappa \left[ x_1^*(\bar{s}) + \delta x_2^*(g(\bar{s}, x_1^*(\bar{s}))) \right] + \delta^2 V_3(g(g(\bar{s}, x_1^*(\bar{s})), x_2^*(g(\bar{s}, x_1^*(\bar{s}))))) \\ &= \ldots \end{align*}Here we show that $V_t$ is differentiable by $s_t$ for all $t=1, 2, \ldots, T+1$
At $t = T+1$, $V_{T+1}(s_{T+1}) = ps_{T+1}$ so this is differentiable by $s_{T+1}$
As an induction hypothesis, assume $V_{t+1}$ is differentiable by $s_t$ ($1 \leq t \leq T$). Since
\begin{align*} V_t(s_t) = -\kappa x_t^*(s_t) + \delta V_{t+1}(g(s_t, x_t^*(s_t))) \end{align*}
and $x_t^*(s_t)$ and $V_{t+1}(s_t)$ are differentiable by $s_t$ and $g(s_t, x)$ is differentiable by both $s_t$ and $x$, $V_t(s_t)$ is also differentiable. The derivative is
\begin{align*} \frac{\partial}{\partial s_t} V_t(s_t) = -\kappa \frac{\partial}{\partial s_t} x_t^*(s) + \delta \frac{\partial}{\partial s_t}V_{t+1}(g(s_t, x_t^*(s_t))) \left\{ g_s(s_t, x_t^*(s_t)) + g_x(s_t, x_t^*(s_t)) \frac{\partial}{\partial s_t} x_t^*(s_t) \right\} \end{align*}
Hence by induction, $V_t$ is differentiable by $s_t$ for all $t$
Next consider the optimality condition of the maximization problem in the Bellman equation: $\underset{x_t \geq 0}{\max} \{ -\kappa x_t + \delta V_{t+1}(g(s_t, x_t)) \}$
Define $\lambda_t(s_t) \equiv V_{t}'(s_t) = \frac{\partial V_t}{\partial s_t}$
Using the chain rule, the FOCs are
\begin{align*} \delta \lambda_{t+1}(g(s_t, x_t^*))g_x(s_t, x_t^*) = \kappa \hspace{1em} t=1, \ldots, T \\ \end{align*}Then the drivatives of the value functions by $s_t$: $\lambda_t(s_t) = \frac{\partial}{\partial s_t} V_t(s_t)$ become
\begin{align*} \lambda_t(s) &= \underbrace{ \bigl\{ -\kappa + \delta \lambda_{t+1}(g(s_t, x_t^*(s_t))) g_x(s_t, x_t^*(s_t)) \bigr\}}_{=0 \text{ by FOC}} \frac{\partial}{\partial s_t} x_t^*(s_t) + \delta \lambda_{t+1}(g(s_t, x_t^*(s_t))) g_s(s_t, x_t^*(s_t)) \\ &= \delta \lambda_{t+1}(g(s_t, x_t^*(s_t))) g_s(s_t, x_t^*(s_t)) \end{align*}for $t=1, \ldots, T$ and
\begin{align*} \lambda_{T+1}(s_{T+1}) = p \end{align*}In summary, the optimal path follows the following equations:
\begin{align*} \begin{cases} \delta \lambda_{t+1}(g(s_t, x_t))g_x(s_t, x_t) = \kappa \hspace{1em} t=1, \ldots, T \\ \lambda_t(s_t) = \delta \lambda_{t+1}(g(s_t, x_t)) g_s(s_t, x_t) \hspace{1em} t=1, \ldots, T \\ \lambda_{T+1}(s_{T+1}) = p \end{cases} \end{align*}How to solve the above equations and get the optimal polity $\{x_t\}_{t=1}^{T}$?
At $t=T$, since $\lambda_{T+1}(s) = p$ regardless of the value of $s$, the 1st equation becomes $$ g_x(s_T, x_T) = \frac{\kappa}{\delta p} $$
So we can get the optimal policy at $t=T$ by solving this equation for $x_T$. Then by the 2nd equation, we get $\lambda_{T}(s_T)$
At $1 \leq t < T$, we can get $x_t$ and $\lambda_t(s_t)$ in the similar way
Finally we set $s_1 = \bar{s}$ and then we get the concrete values of $\lambda_1(s_1^*), x_1^*, \lambda_2(s_2^*), x_2^*, \ldots$
Let
Then
At $t=6$,
At $1 \leq t \leq 5$,
So the optimal feeding policy is
\begin{align*} x_t^*(s_t) = \cfrac{(\delta^{7-t} \alpha^{6-t})^2}{4 \kappa^2} \hspace{1em} t=1, \ldots, 6 \end{align*}Note that in this special case the optimal policy $x_1, \ldots, x_6$ does not depend on the initial weight $\bar{s}$
In addition to the settings of the example above, assume
In order to calculate the value function $V_t(s_t) = \underset{x_t \geq 0}{\max} \{ -\kappa x_t + \delta V_{t+1}(g(s_t, x_t)) \}$ given $s_t$ computationally, we approximate $V_t(s_t)$ as
\begin{align*} V_t(s_t) \simeq \tilde{V_t}(s_t) = \sum_{i=1}^n c_{it} \phi_i(s_t) \end{align*}where $\{ \phi_i \}_{i=1}^n$ are the predetermined basis functions and $\{ c_{it} \}_{i=1}^n$ are coefficients of them
Also we choose $n$ sample points $\xi_1 < \xi_2 < \ldots < \xi_n$ in the state space
The procedure of calculating optimal policy is as follows:
At $t=T$
Using optimizer, solve $V_T(\xi_j) = \underset{x_T \geq 0}{\max} \{ -\kappa x_T + \delta p g(\xi_j, x_T) \}$ for each $\xi_1, \ldots, \xi_n$ and get optimal $V_T(\xi_j)$ and $x_T(\xi_j)$
Solve $n$ simultaneous linear equations
\begin{align*} \sum_{i=1}^n c_{iT} \phi_i(\xi_j) = V_T(\xi_j) \hspace{1em} j=1, \ldots, n \end{align*}
to determine the $n$ coefficients $\{ c_{iT} \}_{i=1}^n$ and get $\tilde{V_T}(s_T)$
At $1 \leq t < T$, given $\tilde{V_{t+1}}(\xi_j)$,
Using optimizer, solve $V_t(\xi_j) = \underset{x_t \geq 0}{\max} \{ -\kappa x_t + \delta \tilde{V_{t+1}}(g(\xi_j, x_t)) \}$ for each $\xi_1, \ldots, \xi_n$ and get optimal $V_t(\xi_j)$ and $x_t(\xi_j)$
Solve $n$ simultaneous linear equations
\begin{align*} \sum_{i=1}^n c_{it} \phi_i(\xi_j) = V_t(\xi_j) \hspace{1em} j=1, \ldots, n \end{align*}
to determine the $n$ coefficients $\{ c_{it} \}_{i=1}^n$ and get $\tilde{V_t}(s_t)$
Finally set $s_1 = \bar{s}$ and then we can get all optimal policies
using QuantEcon
using Optim
using BasisMatrices
using Plots
plotlyjs()
Plotly javascript loaded.
To load again call
init_notebook(true)
Plots.PlotlyJSBackend()
type LiveStock
T::Int # ts_length(sell a livestock at period T+1)
f::Function # cost function for each period
g::Function # state transition function
v_term::Function # value function at period T+1
delta::Float64 # discount factor
basis::Basis # BasisMatrix of value functions
end
# For BasisMatrices, see https://github.com/QuantEcon/BasisMatrices.jl
basis_size = 50
s_min, s_max = 0.4, 2.0
grid = linspace(s_min, s_max, basis_size-2)
basis = Basis(SplineParams(grid, 0, 3)) # cubic spline
1 dimensional Basis on the hypercube formed by (0.4,) × (2.0,). Basis families are Spline
T = 6
delta = 0.9
kappa = 0.4
alpha = 0.9
beta = 0.5
p = 1
f(s, x) = -1 * kappa .* x
g(s, x) = alpha.*s + x.^beta
v_term(s) = p * s
v_term (generic function with 1 method)
live_stock = LiveStock(T, f, g, v_term, delta, basis)
LiveStock(6, f, g, v_term, 0.9, 1 dimensional Basis on the hypercube formed by (0.4,) × (2.0,). Basis families are Spline )
function backward_induction(obj::LiveStock, x_min=0.01, x_max=100)
grid = nodes(obj.basis)[1]
grid_size = length(grid)
VF_coefs = Matrix{Float64}(grid_size, obj.T)
policies = similar(VF_coefs)
Psi = BasisMatrix(obj.basis, Expanded(), grid, 0)
values = Vector{Float64}(grid_size)
# backward induction
for t in obj.T:-1:1
for i in 1:grid_size
if t == obj.T
# minimization -> maximization
obj_func1(x) = -1 * ( obj.f(grid[i], x) + delta * v_term(obj.g(grid[i], x)) )
result = Optim.optimize(obj_func1, x_min, x_max)
else
# minimization -> maximization
obj_func2(x) = -1 * (
obj.f(grid[i], x) +
delta * funeval(VF_coefs[:, t+1], basis, obj.g(grid[i], x))
)
result = Optim.optimize(obj_func2, x_min, x_max)
end
# optimal value
values[i] = -1 * Optim.minimum(result)
# optimal policy
policies[i, t] = Optim.minimizer(result)
end
# calculate coefficients
VF_coefs[:, t] = Psi.vals[1] \ values
end
return VF_coefs, policies
end
backward_induction (generic function with 3 methods)
Let $\bar{s} = 0.6, 1.0, 1.4$
VF_coefs, policies = backward_induction(live_stock)
optimal_values = Vector{Float64}(T+1)
optimal_policy = Vector{Float64}(T)
optimal_livestock_weight = Vector{Float64}(T+1)
initial_s = [0.6, 1.0, 1.4]
for i in 1:length(initial_s)
current_s = initial_s[i]
for t in 1:T
optimal_livestock_weight[t] = current_s
optimal_policy[t] = funeval(policies[:, t], live_stock.basis, current_s)
optimal_values[t] = funeval(VF_coefs[:, t], live_stock.basis, current_s)
current_s = g(current_s, optimal_policy[t])
end
optimal_livestock_weight[T+1] = current_s
optimal_values[T+1] = v_term(current_s)
println("s_init=", initial_s[i])
println(" VF: ", optimal_values)
println(" PF: ", optimal_policy)
println(" We: ", optimal_livestock_weight)
end
s_init=0.6 VF: [1.10697, 1.29836, 1.54685, 1.87759, 2.32835, 2.95458, 3.84596] PF: [0.15387, 0.234523, 0.35745, 0.54481, 0.826924, 1.26697] We: [0.6, 0.932263, 1.32331, 1.78885, 2.34808, 3.02263, 3.84596] s_init=1.0 VF: [1.21995, 1.42388, 1.68632, 2.03256, 2.50053, 3.13976, 4.05219] PF: [0.15387, 0.234523, 0.357449, 0.544786, 0.813225, 1.26804] We: [1.0, 1.29226, 1.64731, 2.08045, 2.6105, 3.25124, 4.05219] s_init=1.4 VF: [1.33292, 1.54941, 1.82579, 2.18753, 2.67192, 3.31609, 4.24875] PF: [0.15387, 0.234522, 0.357449, 0.543014, 0.78211, 1.26948] We: [1.4, 1.65226, 1.97131, 2.37205, 2.87174, 3.46894, 4.24875]
plot feeding policy and livestock weight ($\bar{s}=0.4$)
optimal_values = Vector{Float64}(T+1)
optimal_policy = Vector{Float64}(T)
optimal_livestock_weight = Vector{Float64}(T+1)
initial_s = 0.4
current_s = initial_s
for t in 1:T
optimal_livestock_weight[t] = current_s
optimal_policy[t] = funeval(policies[:, t], live_stock.basis, current_s)
optimal_values[t] = funeval(VF_coefs[:, t], live_stock.basis, current_s)
current_s = g(current_s, optimal_policy[t])
end
optimal_livestock_weight[T+1] = current_s
optimal_values[T+1] = v_term(current_s)
plot(1:7, optimal_livestock_weight, xlims=(0.9, 7.1), ylim=(0, 4), label="(a) Optimal Stock Weight", size=(800, 400))
xlabel!("Month")
ylabel!("Stock Weight")
plot(1:6, optimal_policy, xlims=(0.9, 6.1), ylim=(0, 1.5), label="(b) Optimal Policy", size=(800, 400))
xlabel!("Month")
ylabel!("Feed")
compare to the theoretical policy
thoretical_policy_f(t) = (delta^(7-t) * alpha^(6-t))^2 / 4 / kappa^2
theoretical_policy = [thoretical_policy_f(t) for t in 1:T]
plot(1:6, optimal_policy, xlims=(0.9, 6.1), ylim=(0, 1.5), label="simulated OP", size=(800, 400), lw=3)
plot!(1:6, theoretical_policy, label="theoretical OP", lw=2)
xlabel!("Month")
ylabel!("Feed")