Plotly javascript loaded.
\n", "To load again call
init_notebook(true)\n", " " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "using QuantEcon\n", "using Plots\n", "plotlyjs();" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "price = 1 # Market price of ore\n", "sbar = 100 # Upper bound of ore stock\n", "beta = 0.9 # Discount rate\n", "n = sbar + 1 # Number of states\n", "m = sbar + 1 # Number of actions\n", "\n", "# Cost function\n", "c(s, x) = (x^2) / (1+s);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Product formulation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This approch sets up the reward array R and the transition probability array Q as a 2-dimensional array of shape (n, m) and a 3-simensional array of shape (n, m, n), respectively, where the reward is set to $−∞$ for infeasible state-action pairs (and the transition probability distribution is arbitrary for those pairs)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Reward array:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "R = Matrix{Float64}(n, m)\n", "for j in 1:m\n", " for i in 1:n\n", " if j <= i\n", " R[i, j] = price*(j-1) -c(i-1, j-1)\n", " else\n", " R[i, j] = -Inf\n", " end\n", " end\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(Degenerate) transition probability array:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "Q = zeros(n, m, n)\n", "for j in 1:m\n", " for i in 1:n\n", " if j <= i\n", " Q[i, j, i-j+1] = 1.0\n", " else\n", " Q[i, j, 1] = 1.0\n", " end\n", " end\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set up the dynamic program as a DiscreteDP instance:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ddp = DiscreteDP(R, Q, beta);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve the optimization problem with the solve method, which by defalut uses the policy iteration algorithm:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Solve the dynamic optimization problem (by policy iteration)\n", "res = solve(ddp, PFI);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The number of iterations:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "9" ], "text/plain": [ "9" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "res.num_iter" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The controlled Markov chain is stored in res.mc. To simulate:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0:100" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "res.mc.state_values = 0:100" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "16-element Array{Int64,1}:\n", " 100\n", " 76\n", " 58\n", " 44\n", " 33\n", " 25\n", " 19\n", " 14\n", " 11\n", " 8\n", " 6\n", " 4\n", " 3\n", " 2\n", " 1\n", " 0" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Simulate the controlled Markov chain\n", "nyrs = 15\n", "spath = simulate(res.mc, nyrs+1, init = sbar + 1)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p1 = plot(0:sbar, res.v, xlabel=\"Stock\", ylabel=\"Value\", title=\"Optimal Value Function\")\n", "p2 = plot(0:sbar, res.sigma, xlabel=\"Stock\", ylabel=\"Extraction\", title=\"Optimal Extraction Policy\")\n", "p3 = plot(0:nyrs+1, spath, xlabel=\"Year\", ylabel=\"Stock\", title=\"Optimal State Path\")\n", "plot(p1, p2, p3, legend = false)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## State-action pairs formulation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This approach assigns the rewards and transition probabilities only to feaslbe state-action pairs, setting up R and Q as a 1-dimensional array of length L and a 2-dimensional array of shape (L, n), respectively." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need the arrays of feasible state and action indices:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "5151" ], "text/plain": [ "5151" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Arrays of feasible state and action indices\n", "s_indices = []\n", "a_indices = []\n", "for j in 1:m\n", " for i in 1:n\n", " if i - j >= 0\n", " push!(s_indices, i)\n", " push!(a_indices, j)\n", " end\n", " end\n", "end\n", "\n", "L = length(s_indices)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "5151-element Array{Any,1}:\n", " 1\n", " 2\n", " 3\n", " 4\n", " 5\n", " 6\n", " 7\n", " 8\n", " 9\n", " 10\n", " 11\n", " 12\n", " 13\n", " ⋮\n", " 100\n", " 101\n", " 98\n", " 99\n", " 100\n", " 101\n", " 99\n", " 100\n", " 101\n", " 100\n", " 101\n", " 101" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s_indices" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "5151-element Array{Any,1}:\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " 1\n", " ⋮\n", " 97\n", " 97\n", " 98\n", " 98\n", " 98\n", " 98\n", " 99\n", " 99\n", " 99\n", " 100\n", " 100\n", " 101" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a_indices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Reward vector:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true }, "outputs": [], "source": [ "R_sa = zeros(L)\n", "for (i , (s, x)) in enumerate(zip(s_indices.-1, a_indices.-1))\n", " R_sa[i] = price*x - c(s, x)\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(Degenerate) transition probability array:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Q_sa = Matrix{Float64}(L,n)\n", "for (i, s) in enumerate(s_indices)\n", " x = a_indices[i]\n", " Q_sa[i, s-x+1] = 1\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set up the dynamic program:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true }, "outputs": [], "source": [ "ddp_sa = DiscreteDP(R_sa, Q_sa, beta, s_indices, a_indices);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve the optimization problem with the solve method, which by defalut uses the policy iteration algorithm:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": true }, "outputs": [], "source": [ "res_sa = solve(ddp_sa, PFI);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Number of iterations:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "9" ], "text/plain": [ "9" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "res_sa.num_iter" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Simulate the controlled markov chain:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "16-element Array{Int64,1}:\n", " 100\n", " 76\n", " 58\n", " 44\n", " 33\n", " 25\n", " 19\n", " 14\n", " 11\n", " 8\n", " 6\n", " 4\n", " 3\n", " 2\n", " 1\n", " 0" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "res_sa.mc.state_values = 0:100\n", "nyrs = 15\n", "spath_sa = simulate(res_sa.mc,nyrs+1,init = sbar + 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Draw the graphs:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p1_sa = plot(0:sbar, res_sa.v, xlabel=\"Stock\", ylabel=\"Value\", title=\"Optimal Value Function\")\n", "p2_sa = plot(0:sbar, res_sa.sigma, xlabel=\"Stock\", ylabel=\"Extraction\", title=\"Optimal Extraction Policy\")\n", "p3_sa = plot(0:nyrs+1, spath_sa, xlabel=\"Year\", ylabel=\"Stock\", title=\"Optimal State Path\")\n", "plot(p1_sa, p2_sa, p3_sa, legend=false)" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.6.0", "language": "julia", "name": "julia-0.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6.0" } }, "nbformat": 4, "nbformat_minor": 2 }