import numpy as np
from numba import vectorize, jit
from quantecon.util import tic, toc
import line_profiler
We will define a function to calculate CRRA utility, given a consumption level $c$ and risk-aversion parameter $\gamma$
$$ u(c) = \frac{c^{1-\gamma} - 1}{1 - \gamma} $$** Scalar function **
def CRRA(c, gamma):
return (c ** (1 - gamma) - 1) / (1 - gamma)
CRRA(10, 1.5)
1.367544467966324
Cannot apply to a list
CRRA([10, 5, 6], 1.5)
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-5-711eb6a75c22> in <module>() ----> 1 CRRA([10, 5, 6], 1.5) <ipython-input-4-76cb91a2a3f3> in CRRA(c, gamma) 1 def CRRA(c, gamma): ----> 2 return (c ** (1 - gamma) - 1) / (1 - gamma) 3 4 CRRA(10, 1.5) TypeError: unsupported operand type(s) for ** or pow(): 'list' and 'float'
** NumPy function**
Use numpy to specify $c$ vector
c = np.linspace(5, 25, 1000000)
tic()
CRRA(c, 1.5)
toc()
TOC: Elapsed: 0.043235063552856445 seconds.
0.043235063552856445
** Numba vectorized function - slight speed improvement **
@vectorize
def CRRA_vec(c, gamma):
return (c ** (1 - gamma) - 1) / (1 - gamma)
CRRA_vec(c, 1.5)
tic()
CRRA_vec(c, 1.5)
toc()
TOC: Elapsed: 0.02863907814025879 seconds.
0.02863907814025879
** Numba vectorized function with automatic parallelization - speed doubles **
We can also specify automatic parallelization - this only generates speed improvements with large enough problems
@vectorize('float64(float64, float64)', target='parallel')
def CRRA_parallel(c, gamma):
return (c ** (1 - gamma) - 1) / (1 - gamma)
# Run for compilation
CRRA_parallel(c, 1.5)
tic()
CRRA_parallel(c, 1.5)
toc()
TOC: Elapsed: 0.015202045440673828 seconds.
0.015202045440673828
** Jitted function - no speed improvement over numpy **
@jit
def CRRA_jit(c, gamma):
return (c ** (1 - gamma) - 1) / (1 - gamma)
# Run for compilation
CRRA_jit(c, 1.5)
tic()
CRRA_jit(c, 1.5)
toc()
TOC: Elapsed: 0.03238701820373535 seconds.
0.03238701820373535
Iterative method - slow
grid = np.linspace(-3, 3, 1000)
m = -np.inf
def f(x, y):
return np.cos(x**2 + y**2) / (1 + x**2 + y**2)
tic()
for x in grid:
for y in grid:
z = f(x, y)
if z > m:
m = z
m
toc()
TOC: Elapsed: 2.508265256881714 seconds.
2.508265256881714
Numpy for vectorization - much better
def f(x, y):
return np.cos(x**2 + y**2) / (1 + x**2 + y**2)
x, y = np.meshgrid(grid, grid)
tic()
np.max(f(x, y))
toc()
TOC: Elapsed: 0.048872947692871094 seconds.
0.048872947692871094
Numba for vectorization - even better
import math
@vectorize
def f_vec(x, y):
return math.cos(x**2 + y**2) / (1 + x**2 + y**2)
# Run once for compilation
np.max(f_vec(x, y))
tic()
np.max(f_vec(x, y))
toc()
TOC: Elapsed: 0.022536039352416992 seconds.
0.022536039352416992
Numba for vectorization with automatic parallization - faster even with this small problem (1000 obs)
@vectorize('float64(float64, float64)', target='parallel')
def f_par(x, y):
return math.cos(x**2 + y**2) / (1 + x**2 + y**2)
tic()
np.max(f_par(x, y))
toc()
TOC: Elapsed: 0.014304876327514648 seconds.
0.014304876327514648