added stuff

sandpit/optimal_saving/.ipynb_checkpoints/opt_saving_sandpit-checkpoint.ipynb

@@ -0,0 +1,138 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Optimal Saving\n",
+    "\n",
+    "#### John Stachurski"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Misc imports"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from numba import jit, prange\n",
+    "from quantecon.util import tic, toc"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Pull our class out:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from optsaving import SavingsProblem"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Generate an instance"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sp = SavingsProblem(x_grid_size=200)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Run and time:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Error at iteration 25 is 0.3170445260721948.\n",
+      "Error at iteration 50 is 0.097044286202145.\n",
+      "Error at iteration 75 is 0.03263913361484683.\n",
+      "Error at iteration 100 is 0.011084662423137104.\n",
+      "Error at iteration 125 is 0.0038070466914348344.\n",
+      "Error at iteration 150 is 0.0013257724222484057.\n",
+      "Error at iteration 175 is 0.00046718272051649024.\n",
+      "Error at iteration 200 is 0.00016601430995422106.\n",
+      "\n",
+      "Converged in 213 iterations.\n",
+      "TOC: Elapsed: 0:00:5.52\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "5.521601438522339"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "tic()\n",
+    "v = sp.compute_fixed_point()\n",
+    "toc()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

sandpit/optimal_saving/opt_saving_sandpit.ipynb

@@ -0,0 +1,138 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Optimal Saving\n",
+    "\n",
+    "#### John Stachurski"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Misc imports"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from numba import jit, prange\n",
+    "from quantecon.util import tic, toc"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Pull our class out:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from optsaving import SavingsProblem"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Generate an instance"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sp = SavingsProblem(x_grid_size=200)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Run and time:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Error at iteration 25 is 0.3170445260721948.\n",
+      "Error at iteration 50 is 0.097044286202145.\n",
+      "Error at iteration 75 is 0.03263913361484683.\n",
+      "Error at iteration 100 is 0.011084662423137104.\n",
+      "Error at iteration 125 is 0.0038070466914348344.\n",
+      "Error at iteration 150 is 0.0013257724222484057.\n",
+      "Error at iteration 175 is 0.00046718272051649024.\n",
+      "Error at iteration 200 is 0.00016601430995422106.\n",
+      "\n",
+      "Converged in 213 iterations.\n",
+      "TOC: Elapsed: 0:00:5.52\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "5.521601438522339"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "tic()\n",
+    "v = sp.compute_fixed_point()\n",
+    "toc()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

sandpit/optimal_saving/optsaving.py

@@ -0,0 +1,121 @@
+"""
+
+A simple optimal savings problem.  The Bellman equation is
+
+    v(x, z) = max_x' { u(R x + w z - x') + β E v(x', z')}
+
+where 0 <= x' <= R x + w z and
+
+    E v(x', z') = Σ_{z'} v(x', z') p(z, z')
+
+We take 
+
+    u(c) = c^{1 - γ} / (1 - γ)
+
+and obtain the transition kernel p by discretizing
+
+    log z' = ρ log z + d + σ η
+    
+using Rouwenhorst's method.
+
+"""
+
+import numpy as np
+import quantecon as qe
+from numba import jit, prange
+
+
+@jit(nopython=True)
+def u(c, γ):
+    return (c + 1e-10)**(1 - γ) / (1 - γ)
+
+
+class SavingsProblem:
+
+    def __init__(self, 
+                 β=0.96,
+                 γ=2.5,
+                 ρ=0.9,
+                 d=0.0,
+                 σ=0.1,
+                 r=0.05,
+                 w=1.0,
+                 z_grid_size=25,
+                 x_grid_size=200,
+                 x_grid_max=10):
+
+        self.β, self.γ = β, γ
+        self.R = 1 + r
+        self.w = w
+        self.z_grid_size, self.x_grid_size = z_grid_size, x_grid_size
+
+        mc = qe.rouwenhorst(z_grid_size, d, σ, ρ)
+
+        self.p = mc.P
+        self.z_grid = np.exp(mc.state_values)
+        self.x_grid = np.linspace(0.0, x_grid_max, x_grid_size)
+
+    def pack_parameters(self):
+        return self.β, self.γ, self.R, self.w, self.p, self.x_grid, self.z_grid
+
+    def compute_fixed_point(self, 
+                            tol=1e-4, 
+                            max_iter=1000, 
+                            verbose=True,
+                            print_skip=25): 
+
+        # Set initial condition
+        v_in = np.ones((self.x_grid_size, self.z_grid_size))
+        v_out = np.empty_like(v_in)
+
+        # Set up loop
+        params = self.pack_parameters()
+        i = 0
+        error = tol + 1
+
+        while i < max_iter and error > tol:
+            T(v_in, v_out, params)
+            error = np.max(np.abs(v_in - v_out))
+            i += 1
+            if i % print_skip == 0:
+                print(f"Error at iteration {i} is {error}.")
+            v_in[:] = v_out
+
+        if i == max_iter: 
+            print("Failed to converge!")
+
+        if verbose and i < max_iter:
+            print(f"\nConverged in {i} iterations.")
+
+        return v_out
+
+
+@jit(nopython=True, parallel=True)
+def T(v, v_out, params):
+
+    n, m = v.shape
+    β, γ, R, w, p, x_grid, z_grid = params
+
+    for j in prange(m):
+        z = z_grid[j]
+
+        for i in range(n):
+            x = x_grid[i]
+
+            # Cash in hand at start of period
+            y = R * x + w * z  
+            # A variable to store largest recorded value
+            max_so_far = - np.inf
+
+            # Find largest x_grid index s.t. x' <= y
+            idx = np.searchsorted(x_grid, y)
+
+            # Step through x' with 0 <= x' <= y, find max
+            for k in range(idx):
+                x_next = x_grid[k]
+                val = u(y - x_next, γ) + β * np.sum(v[k, :] * p[j, :])
+                max_so_far = max(val, max_so_far)
+
+            v_out[i, j] = max_so_far
+
+