quant-econ Solutions: Need for Speed

Solutions for http://quant-econ.net/py/need_for_speed.html

Exercise 1

Let's start with some imports



In [1]:

    
import matplotlib.pyplot as plt
import numpy as np
from numba import jit

We let

0 represent "low"
1 represent "high"



In [2]:

    
p, q = 0.1, 0.2  # Prob of leaving low and high state respectively

Here's a pure Python version of the function



In [3]:

    
def compute_series(n):
    x = np.empty(n, dtype=int)
    x[0] = 1  # Start in state 1
    U = np.random.uniform(0, 1, size=n)
    for t in range(1, n):
        current_x = x[t-1]
        if current_x == 0:
            x[t] = U[t] < p
        else:
            x[t] = U[t] > q
    return x

Let's run this code and check that the fraction of time spent in the low state is about 0.666



In [4]:

    
n = 100000
x = compute_series(n)
print(np.mean(x == 0))  # Fraction of time x is in state 0

Now let's time it



In [5]:

    
%timeit compute_series(n)









    



1 loops, best of 3: 198 ms per loop

Next let's implement a Numba version, which is easy



In [6]:

    
compute_series_numba = jit(compute_series)

Let's check we still get the right numbers



In [7]:

    
x = compute_series_numba(n)
print(np.mean(x == 0))

Let's see the time



In [8]:

    
%timeit compute_series_numba(n)









    



1000 loops, best of 3: 1.77 ms per loop

This is a nice speed improvement for one line of code

Now let's implement a Cython version



In [9]:

    
%load_ext cythonmagic



In [10]:

    
%%cython
import numpy as np
from numpy cimport int_t, float_t

def compute_series_cy(int n):
    # == Create NumPy arrays first == #
    x_np = np.empty(n, dtype=int)
    U_np = np.random.uniform(0, 1, size=n)
    # == Now create memoryviews of the arrays == #
    cdef int_t [:] x = x_np
    cdef float_t [:] U = U_np
    # == Other variable declarations == #
    cdef float p = 0.1
    cdef float q = 0.2
    cdef int t
    # == Main loop == #
    x[0] = 1  
    for t in range(1, n):
        current_x = x[t-1]
        if current_x == 0:
            x[t] = U[t] < p
        else:
            x[t] = U[t] > q
    return np.asarray(x)



In [11]:

    
compute_series_cy(10)









    Out[11]:





array([1, 0, 0, 0, 0, 0, 0, 0, 0, 0])



In [12]:

    
x = compute_series_cy(n)
print(np.mean(x == 0))



In [13]:

    
%timeit compute_series_cy(n)









    



100 loops, best of 3: 3.39 ms per loop

The Cython implementation is fast, but not as fast as Numba