{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# DiscreteDP Example: Asset Replacement with Maintenance" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "** Naoyuki Tsuchiya **\n", "\n", "*Department of Economics, University of Tokyo*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From Miranda and Fackler, Applied Computational Economics and Finance, 2002, Section 7.6.3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`Julia` translation of the [Python version](http://nbviewer.jupyter.org/github/QuantEcon/QuantEcon.notebooks/blob/master/ddp_ex_MF_7_6_3_py.ipynb)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "maxage = 5 # Maximum asset age\n", "repcost = 75 # Replacement cost\n", "mancost = 10 # Maintainance cost\n", "beta = 0.9 # Discount factor\n", "m = 3 # Number of actions; 1: keep, 2: service, 3: replace ;" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "25-element Array{Tuple{Int64,Int64},1}:\n", " (1, 0)\n", " (1, 1)\n", " (1, 2)\n", " (1, 3)\n", " (1, 4)\n", " (2, 0)\n", " (2, 1)\n", " (2, 2)\n", " (2, 3)\n", " (2, 4)\n", " (3, 0)\n", " (3, 1)\n", " (3, 2)\n", " (3, 3)\n", " (3, 4)\n", " (4, 0)\n", " (4, 1)\n", " (4, 2)\n", " (4, 3)\n", " (4, 4)\n", " (5, 0)\n", " (5, 1)\n", " (5, 2)\n", " (5, 3)\n", " (5, 4)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "S = vec([(j, i) for i = 0:maxage-1, j = 1:maxage])\n", "S   # State space" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "25" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n = length(S) # Number of states" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "getindex_0 (generic function with 1 method)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# We need a routine to get the index of a age-serv pair\n", "function getindex_0(age, serv, S)\n", " for i in 1:n\n", " if S[i][1] == age && S[i][2] == serv\n", " return i\n", " end\n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "p (generic function with 1 method)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Profit function as a function of the age and the number of service\n", "function p(age, serv)\n", " return (1 - (age - serv)/5) * (50 - 2.5 * age - 2.5 * age^2)\n", "end" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "25×3 Array{Float64,2}:\n", " 36.0 35.0 -25.0\n", " 45.0 44.0 -25.0\n", " 54.0 53.0 -25.0\n", " 63.0 62.0 -25.0\n", " 72.0 71.0 -25.0\n", " 21.0 18.0 -25.0\n", " 28.0 25.0 -25.0\n", " 35.0 32.0 -25.0\n", " 42.0 39.0 -25.0\n", " 49.0 46.0 -25.0\n", " 8.0 2.0 -25.0\n", " 12.0 6.0 -25.0\n", " 16.0 10.0 -25.0\n", " 20.0 14.0 -25.0\n", " 24.0 18.0 -25.0\n", " 0.0 -10.0 -25.0\n", " 0.0 -10.0 -25.0\n", " 0.0 -10.0 -25.0\n", " 0.0 -10.0 -25.0\n", " 0.0 -10.0 -25.0\n", " -Inf -Inf -25.0\n", " -Inf -Inf -25.0\n", " -Inf -Inf -25.0\n", " -Inf -Inf -25.0\n", " -Inf -Inf -25.0" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Reward array\n", "R = Array{Float64}(n, m)\n", "for i in 1:n\n", " R[i, 1] = p(S[i][1], S[i][2])\n", " R[i, 2] = p(S[i][1], S[i][2]+1) - mancost\n", " R[i, 3] = p(0, 0) - repcost\n", "end\n", "\n", "# Infeasible actions\n", "for j in n-maxage+1:n\n", " R[j, 1] = -Inf\n", " R[j, 2] = -Inf\n", "end\n", "R" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# (Degenerate) transition probability array\n", "Q = zeros(n, m, n)\n", "for i in 1:n\n", " Q[i, 1, getindex_0(min(S[i][1]+1, maxage), S[i][2], S)] = 1\n", " Q[i, 2, getindex_0(min(S[i][1]+1, maxage), min(S[i][2]+1, maxage-1), S)] = 1\n", " Q[i, 3, getindex_0(1, 0, S)] = 1\n", "end" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "using QuantEcon" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Create a DiscreteDP\n", "ddp = DiscreteDP(R, Q, beta);" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Solve the dynamic optimization problem (by policy iteration)\n", "res = solve(ddp, PFI);" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "3" ], "text/plain": [ "3" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Number of iterations\n", "res.num_iter" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "25-element Array{Int64,1}:\n", " 2\n", " 2\n", " 2\n", " 2\n", " 1\n", " 1\n", " 2\n", " 2\n", " 2\n", " 1\n", " 3\n", " 1\n", " 1\n", " 1\n", " 1\n", " 3\n", " 3\n", " 3\n", " 3\n", " 3\n", " 3\n", " 3\n", " 3\n", " 3\n", " 3" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Optimal policy\n", "res.sigma" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(1, 0)2(2, 0)1(2, 1)2(3, 0)3(3, 1)1(3, 2)1(4, 0)3(4, 1)3(4, 2)3(4, 3)3(5, 0)3(5, 1)3(5, 2)3(5, 3)3(5, 4)3" ] } ], "source": [ "# Optimal actions for reachable states\n", "for i in 1:n\n", " if S[i][1] > S[i][2]\n", " print(S[i], res.sigma[i])\n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Discrete Markov Chain\n", "stochastic matrix of type Array{Float64,2}:\n", "[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 1.0 0.0 … 0.0 0.0; 1.0 0.0 … 0.0 0.0]" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Simulate the controlled Markov chain\n", "mc = MarkovChain(res.mc.p, S)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "13-element Array{Tuple{Int64,Int64},1}:\n", " (1, 0)\n", " (2, 1)\n", " (3, 2)\n", " (4, 2)\n", " (1, 0)\n", " (2, 1)\n", " (3, 2)\n", " (4, 2)\n", " (1, 0)\n", " (2, 1)\n", " (3, 2)\n", " (4, 2)\n", " (1, 0)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "initial_state_value = getindex_0(1, 0, S)\n", "nyrs = 12\n", "spath = simulate(mc, nyrs+1, init=initial_state_value)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true }, "outputs": [], "source": [ "using Plots" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", " \n", "

Plotly javascript loaded.

\n", "

To load again call

init_notebook(true)

\n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "Plots.PlotlyJSBackend()" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plotlyjs()" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": true }, "outputs": [], "source": [ "spath_1 = Vector{Int}(nyrs)\n", "for i in 1:nyrs\n", " spath_1[i] = spath[i][1]\n", "end\n", "\n", "spath_2 = Vector{Int}(nyrs)\n", "for j in 1:nyrs\n", " spath_2[j] = spath[j][2]\n", "end" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y = spath_1\n", "plot(0:nyrs, y)\n", "title!(\"Optimal State Path: Age of Asset\")\n", "xaxis!(\"Year\")\n", "yaxis!(\"Age of Asset\", 1:0.5:4)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y = spath_2\n", "plot(0:nyrs, y)\n", "title!(\"Optimal State Path: Number of Servicings\")\n", "xaxis!(\"Year\")\n", "yaxis!(\"Number of Servicings\", 0:0.5:2)" ] } ], "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 }