Plotly javascript loaded.
\n", "To load again call
init_notebook(true)\n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "Plots.PlotlyJSBackend()" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using QuantEcon\n", "using Plots\n", "plotlyjs()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "0.9" ], "text/plain": [ "0.9" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The maximum capacity\n", "maxcap = 30\n", "\n", "# The number of states\n", "n = maxcap + 1\n", "\n", "# The number of actions\n", "m = n\n", "\n", "# Constants\n", "a1, b1 = 14, 0.8\n", "a2, b2 = 10, 0.4\n", "\n", "# Reward functions\n", "function F(x::Number)\n", " return a1 * x^b1\n", "end\n", "\n", "function G(x::Number)\n", " return a2 * x^b2\n", "end\n", "\n", "# State transition probabilities\n", "probs = [0.1, 0.2, 0.4, 0.2, 0.1]\n", "supp_size = length(probs)\n", "\n", "# Discount factor\n", "delta = 0.9" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Product formulation" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Reward array\n", "\n", "R = zeros(Float64, n, m)\n", "for s in 1:n\n", " for x in 1:m\n", " \n", " # Because julian index of array begins with 1\n", " xv = x - 1\n", " sv = s - 1\n", " \n", " if x <= s\n", " r = F(xv) + G(sv - xv)\n", " else\n", " r = -Inf\n", " end\n", " \n", " R[s, x] = r\n", " \n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Transition probability array\n", "\n", "Q = zeros(Float64, n, m, n)\n", "for s in 1:n\n", " for x in 1:m\n", " \n", " # Guarding\n", " if x > s\n", " continue\n", " end\n", " \n", " for j in 1:supp_size\n", " Q[s, x, min(s - x + j, n)] += probs[j]\n", " end\n", " \n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "QuantEcon.DPSolveResult{QuantEcon.PFI,Float64}([338.416, 361.988, 377.426, 391.426, 405.236, 417.962, 429.861, 441.171, 452.067, 462.66 … 576.471, 585.14, 593.698, 602.157, 610.511, 618.805, 627.056, 635.21, 643.282, 651.278], [338.416, 361.988, 377.426, 391.426, 405.236, 417.962, 429.861, 441.171, 452.067, 462.66 … 576.471, 585.14, 593.698, 602.157, 610.511, 618.805, 627.056, 635.21, 643.282, 651.278], 4, [1, 1, 1, 2, 2, 2, 2, 2, 2, 2 … 5, 5, 5, 5, 5, 6, 6, 6, 6, 6], Discrete Markov Chain\n", "stochastic matrix of type Array{Float64,2}:\n", "[0.1 0.2 … 0.0 0.0; 0.0 0.1 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.1 0.0])" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solve the dynamic optimization problem by policy iteration\n", "\n", "res = solve(DiscreteDP(R, Q, delta), PFI)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "4" ], "text/plain": [ "4" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Number of iterations\n", "\n", "res.num_iter" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1×31 RowVector{Int64,Array{Int64,1}}:\n", " 1 1 1 2 2 2 2 2 2 2 3 3 3 … 4 4 5 5 5 5 5 6 6 6 6 6" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Optimal policy\n", "\n", "res.sigma'" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1×31 RowVector{Float64,Array{Float64,1}}:\n", " 338.416 361.988 377.426 391.426 … 627.056 635.21 643.282 651.278" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Optimal value function\n", "\n", "res.v'" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1×51 RowVector{Float64,Array{Float64,1}}:\n", " 0.9999 3.004 4.69893 5.89881 … 13.5211 13.513 13.5308 13.5274" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Simulate the controlled Markov chain for num_rep times\n", "# and compute the average\n", "\n", "init = 1\n", "nyrs = 50\n", "ts_length = nyrs + 1\n", "num_rep = 10^4\n", "ave_path = zeros(ts_length)\n", "\n", "for i in 1:num_rep\n", " path = simulate(res.mc, ts_length, init = init)\n", " ave_path = (i / (i + 1)) * ave_path + (1 / (i + 1)) * path\n", "end\n", "\n", "ave_path'" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1×31 RowVector{Float64,Array{Float64,1}}:\n", " 0.0 0.0 1.15115e-7 1.03604e-6 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Stationary distribution of the Markov chain\n", "\n", "stationary_dist = stationary_distributions(res.mc)[1]\n", "\n", "stationary_dist'" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Plot sigma, v, ave_path, stationary_dist\n", "\n", "indices = 0:maxcap\n", "p1 = plot(\n", " scatter(indices, res.sigma .- 1; label=\"Irrigation\"),\n", " title=\"Optimal Irrigation Policy\",\n", " ylabel=\"Irrigation\",\n", " xlabel=\"Water Level\"\n", ")\n", "p2 = plot(\n", " plot(indices, res.v; label=\"Value\"),\n", " title=\"Optimal Value Function\",\n", " ylabel=\"Value\",\n", " xlabel=\"Water Level\"\n", ")\n", "p3 = plot(\n", " plot(ave_path .- 1; label=\"Water Level\"),\n", " title=\"Average Optimal State Path\",\n", " ylabel=\"Water Level\",\n", " xlabel=\"Year\"\n", ")\n", "p4 = plot(\n", " bar(indices, stationary_dist; label=\"Probability\"),\n", " title=\"Stationary Distribution\",\n", " ylabel=\"Probability\",\n", " xlabel=\"Water Level\"\n", ")\n", "\n", "plot(p1, p2, p3, p4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## State-action pairs formulation" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Arrays of state and action indices\n", "\n", "S = 0:n-1\n", "X = 0:m-1\n", "\n", "s_indices = Int64[]\n", "a_indices = Int64[]\n", "S_left = reshape(S, n, 1) .- reshape(S, 1, n)\n", "\n", "for i in 1:n\n", " for j in 1:n\n", " if S_left[i, j] >= 0\n", " push!(s_indices, i)\n", " push!(a_indices, j)\n", " end\n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "496-element Array{Float64,1}:\n", " 0.0 \n", " 10.0 \n", " 14.0 \n", " 13.1951\n", " 24.0 \n", " 24.3754\n", " 15.5185\n", " 27.1951\n", " 34.3754\n", " 33.7151\n", " 17.411 \n", " 29.5185\n", " 37.5705\n", " ⋮ \n", " 173.71 \n", " 178.917 \n", " 184.003 \n", " 188.958 \n", " 193.772 \n", " 198.426 \n", " 202.893 \n", " 207.128 \n", " 211.051 \n", " 214.5 \n", " 217.036 \n", " 212.728 " ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Reward vector\n", "\n", "S_left_o = S_left\n", "S_left = Array{Int64}(length(a_indices))\n", "\n", "for i in 1:length(a_indices)\n", " S_left[i] = S_left_o[s_indices[i], a_indices[i]]\n", "end\n", "\n", "R = F.(X[a_indices]) + G.(S_left)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Transition probability array\n", "\n", "L = length(S_left)\n", "Q = zeros(L, n)\n", "\n", "for i in 1:length(S_left)\n", " s_left = S_left[i]\n", " \n", " for j in 1:supp_size\n", " Q[i, min(s_left + j, n)] += probs[j]\n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "QuantEcon.DPSolveResult{QuantEcon.PFI,Float64}([338.416, 361.988, 377.426, 391.426, 405.236, 417.962, 429.861, 441.171, 452.067, 462.66 … 576.471, 585.14, 593.698, 602.157, 610.511, 618.805, 627.056, 635.21, 643.282, 651.278], [338.416, 361.988, 377.426, 391.426, 405.236, 417.962, 429.861, 441.171, 452.067, 462.66 … 576.471, 585.14, 593.698, 602.157, 610.511, 618.805, 627.056, 635.21, 643.282, 651.278], 4, [1, 1, 1, 2, 2, 2, 2, 2, 2, 2 … 5, 5, 5, 5, 5, 6, 6, 6, 6, 6], Discrete Markov Chain\n", "stochastic matrix of type Array{Float64,2}:\n", "[0.1 0.2 … 0.0 0.0; 0.0 0.1 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.1 0.0])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solve the dynamic optimization problem by policy iteration\n", "\n", "res = solve(DiscreteDP(R, Q, delta, s_indices, a_indices), PFI)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "4" ], "text/plain": [ "4" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Number of iterations\n", "\n", "res.num_iter" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "51-element Array{Float64,1}:\n", " 0.9999 \n", " 3.0158 \n", " 4.68993\n", " 5.88621\n", " 6.93361\n", " 7.94041\n", " 8.83802\n", " 9.59694\n", " 10.2172 \n", " 10.718 \n", " 11.1322 \n", " 11.467 \n", " 11.7598 \n", " ⋮ \n", " 13.4258 \n", " 13.4576 \n", " 13.4336 \n", " 13.4571 \n", " 13.454 \n", " 13.4666 \n", " 13.4685 \n", " 13.4806 \n", " 13.5008 \n", " 13.5062 \n", " 13.5035 \n", " 13.5036 " ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Simulate the controlled Markov chain for num_rep times\n", "# and compute the average\n", "\n", "init = 1\n", "nyrs = 50\n", "ts_length = nyrs + 1\n", "num_rep = 10 ^ 4\n", "ave_path = zeros(ts_length)\n", "\n", "for i in 1:num_rep\n", " path = simulate(res.mc, ts_length, init = init)\n", " ave_path = (i / (i + 1)) * ave_path + (1 / (i + 1)) * path\n", "end\n", "\n", "ave_path" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "31-element Array{Float64,1}:\n", " 0.0 \n", " 0.0 \n", " 1.15115e-7 \n", " 1.03604e-6 \n", " 8.05805e-6 \n", " 5.98598e-5 \n", " 0.000444344\n", " 0.00329805 \n", " 0.0244792 \n", " 0.0608112 \n", " 0.120882 \n", " 0.139769 \n", " 0.15025 \n", " ⋮ \n", " 0.000443927\n", " 5.70763e-5 \n", " 6.34181e-6 \n", " 0.0 \n", " 0.0 \n", " 0.0 \n", " 0.0 \n", " 0.0 \n", " 0.0 \n", " 0.0 \n", " 0.0 \n", " 0.0 " ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Stationary distribution of the Markov chain\n", "\n", "stationary_dist = stationary_distributions(res.mc)[1]" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Plot sigma, v, ave_path, stationary_dist\n", "\n", "indices = 0:maxcap\n", "p1 = plot(\n", " scatter(indices, res.sigma .- 1; label=\"Irrigation\"),\n", " title=\"Optimal Irrigation Policy\",\n", " ylabel=\"Irrigation\",\n", " xlabel=\"Water Level\"\n", ")\n", "p2 = plot(\n", " plot(indices, res.v; label=\"Value\"),\n", " title=\"Optimal Value Function\",\n", " ylabel=\"Value\",\n", " xlabel=\"Water Level\"\n", ")\n", "p3 = plot(\n", " plot(ave_path .- 1; label=\"Water Level\"),\n", " title=\"Average Optimal State Path\",\n", " ylabel=\"Water Level\",\n", " xlabel=\"Year\"\n", ")\n", "p4 = plot(\n", " bar(indices, stationary_dist; label=\"Probability\"),\n", " title=\"Stationary Distribution\",\n", " ylabel=\"Probability\",\n", " xlabel=\"Water Level\"\n", ")\n", "\n", "plot(p1, p2, p3, p4)" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.6.1", "language": "julia", "name": "julia-0.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6.1" } }, "nbformat": 4, "nbformat_minor": 2 }