{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# LQ Approximation with `QuantEcon.jl` and `ContinuousDPs`"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"using QuantEcon\n",
"using ContinuousDPs\n",
"using Plots\n",
"import QuantEcon.ScalarOrArray"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We consider a dynamic maximization problem with\n",
"\n",
"* reward function $f(s, x)$,\n",
"* state transition function $g(s, x)$, and\n",
"* discount rate $\\delta$,\n",
"\n",
"where $s$ and $x$ are the state and the control variables, respectively\n",
"(we follow Miranda-Fackler in notation).\n",
"\n",
"Let $(s^*, x^*)$ denote the steady state state-control pair,\n",
"and write\n",
"$f^* = f(s^*, x^*)$, $f_i^* = f_i(s^*, x^*)$, $f_{ij}^* = f_{ij}(s^*, x^*)$,\n",
"$g^* = g(s^*, x^*)$, and $g_i^* = g_i(s^*, x^*)$ for $i, j = s, x$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First-order expansion of $g$ around $(s^*, x^*)$:\n",
"$$\n",
"\\begin{eqnarray*}\n",
"g\\left(s,x\\right) & \\approx & \\underbrace{g^{*}}_{n\\times1}+\\underbrace{g_{s}^{*}}_{n\\times n}\\left(s-s^{*}\\right)+\\underbrace{g_{x}^{*}}_{n\\times m}\\left(x-x^{*}\\right)\\\\\n",
"\\left[\\begin{array}{c}\n",
"1\\\\\n",
"g\\left(s,x\\right)\n",
"\\end{array}\\right] & \\approx & A\\left(\\begin{array}{c}\n",
"1\\\\\n",
"s\n",
"\\end{array}\\right)+Bx\n",
"\\end{eqnarray*}\n",
"$$\n",
"where\n",
"\n",
"\\begin{eqnarray*}\n",
"A & = & \\left[\\begin{array}{cc}\n",
"1 & 0\\\\\n",
"g^{*}-Dg^{*}z^{*} & g_{s}^{*}\n",
"\\end{array}\\right]\\\\\n",
"B & = & \\left[\\begin{array}{c}\n",
"0\\\\\n",
"g_{x}^{*}\n",
"\\end{array}\\right]\n",
"\\end{eqnarray*}\n",
"\n",
"with $z^* = (s^*, x^*)^{\\mathrm{T}}$ and $Dg^* = (g_s^*, g_x^*)$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Second-order expansion of $f$ around $(s^*, x^*)$:\n",
"\n",
"\\begin{eqnarray*}\n",
"f\\left(s,x,e\\right) & \\approx & \\underbrace{f^{*}}_{1\\times1}+\\underbrace{f_{s}^{*}}_{1\\times n}\\left(s-s^{*}\\right)+\\underbrace{f_{x}^{*}}_{1\\times m}\\left(x-x^{*}\\right)+0.5\\left(s-s^{*}\\right)'\\underbrace{f_{ss}^{*}}_{n\\times n}\\left(s-s^{*}\\right)\\\\\n",
" & & +\\left(s-s^{*}\\right)'\\underbrace{f_{sx}^{*}}_{n\\times m}\\left(x-x^{*}\\right)+0.5\\left(x-x^{*}\\right)'\\underbrace{f_{xx}^{*}}_{m\\times m}\\left(x-x^{*}\\right)\\\\\n",
" & = & -\\left(\\left(\\begin{array}{cc}\n",
"1 & s'\\end{array}\\right)R\\left(\\begin{array}{c}\n",
"1\\\\\n",
"s\n",
"\\end{array}\\right)+2x'N\\left(\\begin{array}{c}\n",
"1\\\\\n",
"s\n",
"\\end{array}\\right)+x'Qx\\right)\n",
"\\end{eqnarray*}\n",
"\n",
"where\n",
"\n",
"\\begin{eqnarray*}\n",
"R & = & -\\left[\\begin{array}{cc}\n",
"f^{*}-\\left(Df^{*}\\right)'z^{*}+0.5\\left(z^{*}\\right)'D^{2}f^{*}z^{*} & 0.5\\left(f_{s}^{*}-\\left(\\left(s^{*}\\right)'f_{ss}^{*}+\\left(f_{sx}^{*}x^{*}\\right)'\\right)\\right)\\\\\n",
"0.5\\left(\\left(f_{s}^{*}\\right)'-\\left(f_{ss}^{*}s^{*}+f_{sx}^{*}x^{*}\\right)\\right) & 0.5f_{ss}^{*}\n",
"\\end{array}\\right]\\\\\n",
"N & = & -0.5\\left[\\begin{array}{cc}\n",
"\\left(\\left(f_{x}^{*}\\right)'-\\left(\\left(f_{sx}^{*}\\right)'s^{*}+f_{xx}^{*}x^{*}\\right)\\right) & \\left(f_{sx}^{*}\\right)'\\end{array}\\right]\\\\\n",
"Q & = & -\\begin{array}{c}\n",
"\\left[0.5f_{xx}^{*}\\right]\\end{array}\n",
"\\end{eqnarray*}\n",
"\n",
"with $Df^{*}=\\left[\\begin{array}{cc}\n",
"f_{s}^{*} & f_{x}^{*}\\end{array}\\right]^\\mathrm{T}$ and $D^{2}f^{*}=\\left[\\begin{array}{cc}\n",
"f_{ss}^{*} & f_{sx}^{*}\\\\\n",
"f_{xs}^{*} & f_{xx}^{*}\n",
"\\end{array}\\right]$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Optimal Economic Growth"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We consider the following optimal growth model from Miranda and Fackler, Section 9.7.1:\n",
"\n",
"* $f(s, x) = \\dfrac{(s - x)^{1-\\alpha}}{1-\\alpha}$,\n",
"* $g(s, x) = \\gamma + x^{\\beta}$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"alpha = 0.2\n",
"bet = 0.5\n",
"gamm = 0.9\n",
"discount = 0.9;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Function definitions:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"f(s, x) = (s - x)^(1 - alpha) / (1 - alpha)\n",
"f_s(s, x) = (s - x)^(-alpha)\n",
"f_x(s, x) = -f_s(s, x)\n",
"f_ss(s, x) = -alpha * (s - x)^(-alpha - 1)\n",
"f_sx(s, x) = -f_ss(s, x)\n",
"f_xx(s, x) = f_ss(s, x)\n",
"\n",
"g(s, x) = gamm * x + x^bet\n",
"g_s(s, x) = 0\n",
"g_x(s, x) = gamm + bet * x^(bet - 1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Steady state:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(7.416897506925212, 5.6094182825484795)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x_star = ((discount * bet) / (1 - discount * gamm))^(1 / (1 - bet))\n",
"s_star = gamm * x_star + x_star^bet\n",
"s_star, x_star"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`(s_star, x_star)` satisfies the Euler equations:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.1102230246251565e-16"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f_x(s_star, x_star) + discount * f_s(g(s_star, x_star), x_star) * g_x(s_star, x_star)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Construct $f^*$, $D f^*$, $D^2 f^*$, $g^*$, and $D g^*$:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"f_star = f(s_star, x_star)\n",
"Df_star = [f_s(s_star, x_star), f_x(s_star, x_star)]\n",
"DDf_star = [f_ss(s_star, x_star) f_sx(s_star, x_star);\n",
" f_sx(s_star, x_star) f_xx(s_star, x_star)]\n",
"g_star = g(s_star, x_star)\n",
"Dg_star = [g_s(s_star, x_star), g_x(s_star, x_star)];"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### LQ Approximation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Generate an LQ instance that approximates our dynamic optimization problem:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"LQ([0.04914869864113181], [-0.24085180775732695 -0.5330115101939921; -0.5330115101939921 0.04914869864113181], [1.0 0.0; 1.1842105263157894 0.0], [0.0; 1.1111111111111112], [0.0; 0.0], [0.5330115101939921 -0.04914869864113181], 0.9, nothing, [NaN NaN; NaN NaN], [NaN NaN; NaN NaN], 0.0, [0.0 0.0])"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lq = approx_lq(s_star, x_star, f_star, Df_star, DDf_star, g_star, Dg_star', discount)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solution by `stationary_values(::LQ)`"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Solve the LQ problem:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"([-13.211795608258779 -0.4806293445369966; -0.4806293445369966 0.004914869864113194], [1.0657894736842055 -0.8999999999999997], 0.0)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P, F, d = stationary_values(lq)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The optimal value function (of the LQ minimization problem):"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"V (generic function with 1 method)"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"V(s) = [1, s]' * P * [1, s] + d"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The value at $s^*$:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"-20.070983979777214"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"V(s_star)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-20.070983979777242"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"-f_star / (1 - lq.bet)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The optimal policy function:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"X (generic function with 1 method)"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X(s) = - (F * [1, s])[1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The optimal choice at $s^*$:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5.609418282548483"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X(s_star)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"5.6094182825484795"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x_star"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s_min, s_max = 5, 10\n",
"ss = range(s_min, stop=s_max, length=50)\n",
"title = \"Optimal Investment Policy\"\n",
"xlabel = \"Wealth\"\n",
"ylabel = \"Investment (% of Wealth)\"\n",
"plot(ss, X.(ss)./ss, xlims=(s_min, s_max), ylims=(0.65, 0.9),\n",
" title=title, xlabel=xlabel, ylabel=ylabel, label=\"L-Q\")\n",
"plot!([s_star], [x_star/s_star], m=(7,:star8), label=\"\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Renewable Resource Management"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider the renewable resource management model from Miranda and Fackler, Section 9.7.2:\n",
"\n",
"* $f(s, x) = \\dfrac{(s - x)^{1-\\gamma}}{1-\\gamma} - \\kappa (s - x)$,\n",
"* $g(s, x) = \\alpha x - 0.5 \\beta x^2$."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"alpha = 4.0\n",
"bet = 1.0\n",
"gamm = 0.5\n",
"kappa = 0.2\n",
"discount = 0.9;"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"f(s, x) = (s - x)^(1 - gamm) / (1 - gamm) - kappa * (s - x)\n",
"f_s(s, x) = (s - x)^(-gamm) - kappa\n",
"f_x(s, x) = -f_s(s, x)\n",
"f_ss(s, x) = -gamm * (s - x)^(-gamm - 1)\n",
"f_sx(s, x) = -f_ss(s, x)\n",
"f_xx(s, x) = f_ss(s, x)\n",
"\n",
"g(s, x) = alpha * x - 0.5 * bet * x^2\n",
"g_s(s, x) = 0\n",
"g_x(s, x) = alpha - bet * x;"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(7.382716049382716, 2.888888888888889)"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x_star = (discount * alpha - 1) / (discount * bet)\n",
"s_star = (alpha^2 - 1/discount^2) / (2 * bet)\n",
"s_star, x_star"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5.551115123125783e-17"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f_x(s_star, x_star) + discount * f_s(g(s_star, x_star), x_star) * g_x(s_star, x_star)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"f_star = f(s_star, x_star)\n",
"Df_star = [f_s(s_star, x_star), f_x(s_star, x_star)]\n",
"DDf_star = [f_ss(s_star, x_star) f_sx(s_star, x_star);\n",
" f_sx(s_star, x_star) f_xx(s_star, x_star)]\n",
"g_star = g(s_star, x_star)\n",
"Dg_star = [g_s(s_star, x_star) g_x(s_star, x_star)];"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"LQ([0.0262431197105178], [-1.589898669028243 -0.2537961323936474; -0.2537961323936474 0.0262431197105178], [1.0 0.0; 4.172839506172839 0.0], [0.0; 1.1111111111111112], [0.0; 0.0], [0.2537961323936474 -0.0262431197105178], 0.9, nothing, [NaN NaN; NaN NaN], [NaN NaN; NaN NaN], 0.0, [0.0 0.0])"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lq = approx_lq(s_star, x_star, f_star, Df_star, DDf_star, g_star, Dg_star, discount)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"([-31.26051474783933 -0.15523863836970273; -0.15523863836970286 0.0026243119710518074], [3.755555555555553 -0.899999999999999], 0.0)"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P, F, d = stationary_values(lq)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"V (generic function with 1 method)"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"V(s) = [1, s]' * P * [1, s] + d"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-33.40964351976545"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"V(s_star)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-33.409643519765496"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"-f_star / (1 - lq.bet)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"X (generic function with 1 method)"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X(s) = - (F * [1, s])[1]"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.888888888888884"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X(s_star)"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.888888888888889"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x_star"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s_min, s_max = 6, 9\n",
"ss = range(s_min, stop=s_max, length=50)\n",
"harvest = ss - X.(ss)\n",
"h_star = s_star - x_star\n",
"title = \"Optimal Harvest Policy\"\n",
"xlabel = \"Available Stock\"\n",
"ylabel = \"Harvest (% of Stock)\"\n",
"plot(ss, harvest./ss, xlims=(s_min, s_max), ylims=(0.5, 0.75),\n",
" title=title, xlabel=xlabel, ylabel=ylabel, label=\"L-Q\")\n",
"plot!([s_star], [h_star/s_star], m=(7,:star8), label=\"\")"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"shadow_price(s) = -2 * (P * [1, s])[2]\n",
"title = \"Shadow Price Function\"\n",
"ylabel = \"Price\"\n",
"plot(ss, shadow_price.(ss), xlims=(s_min, s_max), ylims=(0.2, 0.4),\n",
" title=title, xlabel=xlabel, ylabel=ylabel, label=\"L-Q\")\n",
"plot!([s_star], [shadow_price(s_star)], m=(7,:star8), label=\"\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.5.3",
"language": "julia",
"name": "julia-1.5"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.5.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}