<br>

# Python for Finance

**Analyze Big Financial Data**

O'Reilly (2014)

Yves Hilpisch



**Buy the book ** |
<a href='http://shop.oreilly.com/product/0636920032441.do' target='_blank'>O'Reilly</a> |
<a href='http://www.amazon.com/Yves-Hilpisch/e/B00JCYHHJM' target='_blank'>Amazon</a>

**All book codes & IPYNBs** |
<a href="http://oreilly.quant-platform.com">http://oreilly.quant-platform.com</a>

**The Python Quants GmbH** | <a href='http://pythonquants.com' target='_blank'>www.pythonquants.com</a>

**Contact us** | <a href='mailto:analytics@pythonquants.com'>analytics@pythonquants.com</a>

# Mathematical Tools

## Approximation


```python
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
```


```python
def f(x):
    return np.sin(x) + 0.5 * x
```


```python
x = np.linspace(-2 * np.pi, 2 * np.pi, 50)
```


```python
plt.plot(x, f(x), 'b')
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')
# tag: sin_plot
# title: Example function plot
# size: 60
```

### Regression

#### Monomials as Basis Functions


```python
reg = np.polyfit(x, f(x), deg=1)
ry = np.polyval(reg, x)
```


```python
plt.plot(x, f(x), 'b', label='f(x)')
plt.plot(x, ry, 'r.', label='regression')
plt.legend(loc=0)
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')
# tag: sin_plot_reg_1
# title: Example function and linear regression
# size: 60
```


```python
reg = np.polyfit(x, f(x), deg=5)
ry = np.polyval(reg, x)
```


```python
plt.plot(x, f(x), 'b', label='f(x)')
plt.plot(x, ry, 'r.', label='regression')
plt.legend(loc=0)
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')
# tag: sin_plot_reg_2
# title: Regression with monomials up to order 5
# size: 60
```


```python
reg = np.polyfit(x, f(x), 7)
ry = np.polyval(reg, x)
```


```python
plt.plot(x, f(x), 'b', label='f(x)')
plt.plot(x, ry, 'r.', label='regression')
plt.legend(loc=0)
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')
# tag: sin_plot_reg_3
# title: Regression with monomials up to order 7
# size: 60
```


```python
np.allclose(f(x), ry)
```




    False




```python
np.sum((f(x) - ry) ** 2) / len(x)
```




    0.0017769134759517413



#### Individual Basis Functions


```python
matrix = np.zeros((3 + 1, len(x)))
matrix[3, :] = x ** 3
matrix[2, :] = x ** 2
matrix[1, :] = x
matrix[0, :] = 1
```


```python
reg = np.linalg.lstsq(matrix.T, f(x))[0]
```


```python
reg
```




    array([  1.13968447e-14,   5.62777448e-01,  -8.88178420e-16,
            -5.43553615e-03])




```python
ry = np.dot(reg, matrix)
```


```python
plt.plot(x, f(x), 'b', label='f(x)')
plt.plot(x, ry, 'r.', label='regression')
plt.legend(loc=0)
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')
# tag: sin_plot_reg_4
# title: Regression via least-squares function
# size: 60
```


```python
matrix[3, :] = np.sin(x)
reg = np.linalg.lstsq(matrix.T, f(x))[0]
ry = np.dot(reg, matrix)
```


```python
plt.plot(x, f(x), 'b', label='f(x)')
plt.plot(x, ry, 'r.', label='regression')
plt.legend(loc=0)
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')
# tag: sin_plot_reg_5
# title: Regression using individual functions
# size: 60
```


```python
np.allclose(f(x), ry)
```




    True




```python
np.sum((f(x) - ry) ** 2) / len(x)
```




    2.2749084503102031e-31




```python
reg
```




    array([  1.55428020e-16,   5.00000000e-01,   0.00000000e+00,
             1.00000000e+00])



#### Noisy Data


```python
xn = np.linspace(-2 * np.pi, 2 * np.pi, 50)
xn = xn + 0.15 * np.random.standard_normal(len(xn))
yn = f(xn) + 0.25 * np.random.standard_normal(len(xn))
```


```python
reg = np.polyfit(xn, yn, 7)
ry = np.polyval(reg, xn)
```


```python
plt.plot(xn, yn, 'b^', label='f(x)')
plt.plot(xn, ry, 'ro', label='regression')
plt.legend(loc=0)
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')
# tag: sin_plot_reg_6
# title: Regression with noisy data
# size: 60
```

#### Unsorted Data


```python
xu = np.random.rand(50) * 4 * np.pi - 2 * np.pi
yu = f(xu)
```


```python
print xu[:10].round(2)
print yu[:10].round(2)
```

    [ 5.22  4.42  2.61  5.37 -4.4  -4.1  -2.64  1.84 -5.12  5.09]
    [ 1.74  1.25  1.81  1.89 -1.25 -1.23 -1.8   1.88 -1.64  1.61]



```python
reg = np.polyfit(xu, yu, 5)
ry = np.polyval(reg, xu)
```


```python
plt.plot(xu, yu, 'b^', label='f(x)')
plt.plot(xu, ry, 'ro', label='regression')
plt.legend(loc=0)
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')
# tag: sin_plot_reg_7
# title: Regression with unsorted data
# size: 60
```

#### Multiple Dimensions


```python
def fm((x, y)):
    return np.sin(x) + 0.25 * x + np.sqrt(y) + 0.05 * y ** 2
```


```python
x = np.linspace(0, 10, 20)
y = np.linspace(0, 10, 20)
X, Y = np.meshgrid(x, y)
  # generates 2-d grids out of the 1-d arrays
Z = fm((X, Y))
x = X.flatten()
y = Y.flatten()
  # yields 1-d arrays from the 2-d grids
```


```python
from mpl_toolkits.mplot3d import Axes3D
import matplotlib as mpl

fig = plt.figure(figsize=(9, 6))
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap=mpl.cm.coolwarm,
        linewidth=0.5, antialiased=True)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x, y)')
fig.colorbar(surf, shrink=0.5, aspect=5)
# tag: sin_plot_3d_1
# title: Function with two parameters
# size: 60
```


```python
matrix = np.zeros((len(x), 6 + 1))
matrix[:, 6] = np.sqrt(y)
matrix[:, 5] = np.sin(x)
matrix[:, 4] = y ** 2
matrix[:, 3] = x ** 2
matrix[:, 2] = y
matrix[:, 1] = x
matrix[:, 0] = 1
```


```python
import statsmodels.api as sm
```


```python
model = sm.OLS(fm((x, y)), matrix).fit()
```


```python
model.rsquared
```




    1.0




```python
a = model.params
a
```




    array([  7.14706072e-15,   2.50000000e-01,  -2.22044605e-16,
            -1.02348685e-16,   5.00000000e-02,   1.00000000e+00,
             1.00000000e+00])




```python
def reg_func(a, (x, y)):
    f6 = a[6] * np.sqrt(y)
    f5 = a[5] * np.sin(x)
    f4 = a[4] * y ** 2
    f3 = a[3] * x ** 2
    f2 = a[2] * y
    f1 = a[1] * x
    f0 = a[0] * 1
    return (f6 + f5 + f4 + f3 +
            f2 + f1 + f0)
```


```python
RZ = reg_func(a, (X, Y))
```


```python
fig = plt.figure(figsize=(9, 6))
ax = fig.gca(projection='3d')
surf1 = ax.plot_surface(X, Y, Z, rstride=2, cstride=2,
            cmap=mpl.cm.coolwarm, linewidth=0.5,
            antialiased=True)
surf2 = ax.plot_wireframe(X, Y, RZ, rstride=2, cstride=2,
                          label='regression')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x, y)')
ax.legend()
fig.colorbar(surf, shrink=0.5, aspect=5)
# tag: sin_plot_3d_2
# title: Higher dimension regression
# size: 60
```

### Interpolation


```python
import scipy.interpolate as spi
```


```python
x = np.linspace(-2 * np.pi, 2 * np.pi, 25)
```


```python
def f(x):
    return np.sin(x) + 0.5 * x
```


```python
ipo = spi.splrep(x, f(x), k=1)
```


```python
iy = spi.splev(x, ipo)
```


```python
plt.plot(x, f(x), 'b', label='f(x)')
plt.plot(x, iy, 'r.', label='interpolation')
plt.legend(loc=0)
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')
# tag: sin_plot_ipo_1
# title: Example plot with linear interpolation
# size: 60
```


```python
np.allclose(f(x), iy)
```




    True




```python
xd = np.linspace(1.0, 3.0, 50)
iyd = spi.splev(xd, ipo)
```


```python
plt.plot(xd, f(xd), 'b', label='f(x)')
plt.plot(xd, iyd, 'r.', label='interpolation')
plt.legend(loc=0)
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')
# tag: sin_plot_ipo_2
# title: Example plot (detail) with linear interpolation
# size: 60
```


```python
ipo = spi.splrep(x, f(x), k=3)
iyd = spi.splev(xd, ipo)
```


```python
plt.plot(xd, f(xd), 'b', label='f(x)')
plt.plot(xd, iyd, 'r.', label='interpolation')
plt.legend(loc=0)
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')
# tag: sin_plot_ipo_3
# title: Example plot (detail) with cubic splines interpolation
# size: 60
```


```python
np.allclose(f(xd), iyd)
```




    False




```python
np.sum((f(xd) - iyd) ** 2) / len(xd)
```




    1.1349319851436252e-08



## Convex Optimization


```python
def fm((x, y)):
    return (np.sin(x) + 0.05 * x ** 2
          + np.sin(y) + 0.05 * y ** 2)
```


```python
x = np.linspace(-10, 10, 50)
y = np.linspace(-10, 10, 50)
X, Y = np.meshgrid(x, y)
Z = fm((X, Y))
```


```python
fig = plt.figure(figsize=(9, 6))
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap=mpl.cm.coolwarm,
        linewidth=0.5, antialiased=True)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x, y)')
fig.colorbar(surf, shrink=0.5, aspect=5)
# tag: opt_plot_3d
# title: Function to minimize with two parameters
# size: 60
```


```python
import scipy.optimize as spo
```

### Global Optimization


```python
def fo((x, y)):
    z = np.sin(x) + 0.05 * x ** 2 + np.sin(y) + 0.05 * y ** 2
    if output == True:
        print '%8.4f %8.4f %8.4f' % (x, y, z)
    return z
```


```python
output = True
spo.brute(fo, ((-10, 10.1, 5), (-10, 10.1, 5)), finish=None)
```

    -10.0000 -10.0000  11.0880
    -10.0000 -10.0000  11.0880
    -10.0000  -5.0000   7.7529
    -10.0000   0.0000   5.5440
    -10.0000   5.0000   5.8351
    -10.0000  10.0000  10.0000
     -5.0000 -10.0000   7.7529
     -5.0000  -5.0000   4.4178
     -5.0000   0.0000   2.2089
     -5.0000   5.0000   2.5000
     -5.0000  10.0000   6.6649
      0.0000 -10.0000   5.5440
      0.0000  -5.0000   2.2089
      0.0000   0.0000   0.0000
      0.0000   5.0000   0.2911
      0.0000  10.0000   4.4560
      5.0000 -10.0000   5.8351
      5.0000  -5.0000   2.5000
      5.0000   0.0000   0.2911
      5.0000   5.0000   0.5822
      5.0000  10.0000   4.7471
     10.0000 -10.0000  10.0000
     10.0000  -5.0000   6.6649
     10.0000   0.0000   4.4560
     10.0000   5.0000   4.7471
     10.0000  10.0000   8.9120





    array([ 0.,  0.])




```python
output = False
opt1 = spo.brute(fo, ((-10, 10.1, 0.1), (-10, 10.1, 0.1)), finish=None)
opt1
```




    array([-1.4, -1.4])




```python
fm(opt1)
```




    -1.7748994599769203



### Local Optimization


```python
output = True
opt2 = spo.fmin(fo, opt1, xtol=0.001, ftol=0.001, maxiter=15, maxfun=20)
opt2
```

     -1.4000  -1.4000  -1.7749
     -1.4700  -1.4000  -1.7743
     -1.4000  -1.4700  -1.7743
     -1.3300  -1.4700  -1.7696
     -1.4350  -1.4175  -1.7756
     -1.4350  -1.3475  -1.7722
     -1.4088  -1.4394  -1.7755
     -1.4438  -1.4569  -1.7751
     -1.4328  -1.4427  -1.7756
     -1.4591  -1.4208  -1.7752
     -1.4213  -1.4347  -1.7757
     -1.4235  -1.4096  -1.7755
     -1.4305  -1.4344  -1.7757
     -1.4168  -1.4516  -1.7753
     -1.4305  -1.4260  -1.7757
     -1.4396  -1.4257  -1.7756
     -1.4259  -1.4325  -1.7757
     -1.4259  -1.4241  -1.7757
     -1.4304  -1.4177  -1.7757
     -1.4270  -1.4288  -1.7757
    Warning: Maximum number of function evaluations has been exceeded.





    array([-1.42702972, -1.42876755])




```python
fm(opt2)
```




    -1.7757246992239009




```python
output = False
spo.fmin(fo, (2.0, 2.0), maxiter=250)
```

    Optimization terminated successfully.
             Current function value: 0.015826
             Iterations: 46
             Function evaluations: 86





    array([ 4.2710728 ,  4.27106945])



### Constrained Optimization


```python
# function to be minimized
from math import sqrt
def Eu((s, b)):
    return -(0.5 * sqrt(s * 15 + b * 5) + 0.5 * sqrt(s * 5 + b * 12))

# constraints
cons = ({'type': 'ineq', 'fun': lambda (s, b):  100 - s * 10 - b * 10})
  # budget constraint
bnds = ((0, 1000), (0, 1000))  # uppper bounds large enough
```


```python
result = spo.minimize(Eu, [5, 5], method='SLSQP',
                       bounds=bnds, constraints=cons)
```


```python
result
```




      status: 0
     success: True
        njev: 5
        nfev: 21
         fun: -9.700883611487832
           x: array([ 8.02547122,  1.97452878])
     message: 'Optimization terminated successfully.'
         jac: array([-0.48508096, -0.48489535,  0.        ])
         nit: 5




```python
result['x']
```




    array([ 8.02547122,  1.97452878])




```python
-result['fun']
```




    9.700883611487832




```python
np.dot(result['x'], [10, 10])
```




    99.999999999999986



## Integration


```python
import scipy.integrate as sci
```


```python
def f(x):
    return np.sin(x) + 0.5 * x
```


```python
a = 0.5  # left integral limit
b = 9.5  # right integral limit
x = np.linspace(0, 10)
y = f(x)
```


```python
from matplotlib.patches import Polygon

fig, ax = plt.subplots(figsize=(7, 5))
plt.plot(x, y, 'b', linewidth=2)
plt.ylim(ymin=0)

# area under the function
# between lower and upper limit
Ix = np.linspace(a, b)
Iy = f(Ix)
verts = [(a, 0)] + list(zip(Ix, Iy)) + [(b, 0)]
poly = Polygon(verts, facecolor='0.7', edgecolor='0.5')
ax.add_patch(poly)

# labels
plt.text(0.75 * (a + b), 1.5, r"$\int_a^b f(x)dx$",
         horizontalalignment='center', fontsize=20)

plt.figtext(0.9, 0.075, '$x$')
plt.figtext(0.075, 0.9, '$f(x)$')

ax.set_xticks((a, b))
ax.set_xticklabels(('$a$', '$b$'))
ax.set_yticks([f(a), f(b)])
# tag: sin_integral
# title: Example function with integral area
# size: 50
```

### Numerical Integration


```python
sci.fixed_quad(f, a, b)[0]
```




    24.366995967084602




```python
sci.quad(f, a, b)[0]
```




    24.374754718086752




```python
sci.romberg(f, a, b)
```




    24.374754718086713




```python
xi = np.linspace(0.5, 9.5, 25)
```


```python
sci.trapz(f(xi), xi)
```




    24.352733271544516




```python
sci.simps(f(xi), xi)
```




    24.374964184550748



### Integration by Simulation


```python
for i in range(1, 20):
    np.random.seed(1000)
    x = np.random.random(i * 10) * (b - a) + a
    print np.sum(f(x)) / len(x) * (b - a)
```

    24.8047622793
    26.5229188983
    26.2655475192
    26.0277033994
    24.9995418144
    23.8818101416
    23.5279122748
    23.507857659
    23.6723674607
    23.6794104161
    24.4244017079
    24.2390053468
    24.115396925
    24.4241919876
    23.9249330805
    24.1948421203
    24.1173483782
    24.1006909297
    23.7690510985


## Symbolic Computation


```python
import sympy as sy
```

### Basics


```python
x = sy.Symbol('x')
y = sy.Symbol('y')
```


```python
type(x)
```




    sympy.core.symbol.Symbol




```python
sy.sqrt(x)
```




    sqrt(x)




```python
3 + sy.sqrt(x) - 4 ** 2
```




    sqrt(x) - 13




```python
f = x ** 2 + 3 + 0.5 * x ** 2 + 3 / 2
```


```python
sy.simplify(f)
```




    1.5*x**2 + 4




```python
sy.init_printing(pretty_print=False, use_unicode=False)
```


```python
print sy.pretty(f)
```

         2    
    1.5*x  + 4



```python
print sy.pretty(sy.sqrt(x) + 0.5)
```

      ___      
    \/ x  + 0.5



```python
pi_str = str(sy.N(sy.pi, 400000))
pi_str[:40]
```




    '3.14159265358979323846264338327950288419'




```python
pi_str[-40:]
```




    '8245672736856312185020980470362464176198'




```python
pi_str.find('111272')
```




    366713



### Equations


```python
sy.solve(x ** 2 - 1)
```




    [-1, 1]




```python
sy.solve(x ** 2 - 1 - 3)
```




    [-2, 2]




```python
sy.solve(x ** 3 + 0.5 * x ** 2 - 1)
```




    [0.858094329496553, -0.679047164748276 - 0.839206763026694*I, -0.679047164748276 + 0.839206763026694*I]




```python
sy.solve(x ** 2 + y ** 2)
```




    [{x: -I*y}, {x: I*y}]



### Integration


```python
a, b = sy.symbols('a b')
```


```python
print sy.pretty(sy.Integral(sy.sin(x) + 0.5 * x, (x, a, b)))
```

      b                    
      /                    
     |                     
     |  (0.5*x + sin(x)) dx
     |                     
    /                      
    a                      



```python
int_func = sy.integrate(sy.sin(x) + 0.5 * x, x)
```


```python
print sy.pretty(int_func)
```

          2         
    0.25*x  - cos(x)



```python
Fb = int_func.subs(x, 9.5).evalf()
Fa = int_func.subs(x, 0.5).evalf()
```


```python
Fb - Fa  # exact value of integral
```




    24.3747547180867




```python
int_func_limits = sy.integrate(sy.sin(x) + 0.5 * x, (x, a, b))
print sy.pretty(int_func_limits)
```

               2         2                  
    - - -0.25*a  + 0.25*b  + cos(a) - cos(b)



```python
int_func_limits.subs({a : 0.5, b : 9.5}).evalf()
```




    24.3747547180868




```python
sy.integrate(sy.sin(x) + 0.5 * x, (x, 0.5, 9.5))
```




    24.3747547180868



### Differentiation


```python
int_func.diff()
```




    0.5*x + sin(x)




```python
f = (sy.sin(x) + 0.05 * x ** 2
   + sy.sin(y) + 0.05 * y ** 2)
```


```python
del_x = sy.diff(f, x)
del_x
```




    0.1*x + cos(x)




```python
del_y = sy.diff(f, y)
del_y
```




    0.1*y + cos(y)




```python
xo = sy.nsolve(del_x, -1.5)
xo
```




    mpf('-1.4275517787645941')




```python
yo = sy.nsolve(del_y, -1.5)
yo
```




    mpf('-1.4275517787645941')




```python
f.subs({x : xo, y : yo}).evalf() 
  # global minimum
```




    -1.77572565314742




```python
xo = sy.nsolve(del_x, 1.5)
xo
```




    mpf('1.7463292822528528')




```python
yo = sy.nsolve(del_y, 1.5)
yo
```




    mpf('1.7463292822528528')




```python
f.subs({x : xo, y : yo}).evalf()
  # local minimum
```




    2.27423381055640



## Conclusions

## Further Reading

<br>

<a href="http://www.pythonquants.com" target="_blank">www.pythonquants.com</a> | <a href="http://twitter.com/dyjh" target="_blank">@dyjh</a>

<a href="mailto:analytics@pythonquants.com">analytics@pythonquants.com</a>

**Python Quant Platform** |
<a href="http://oreilly.quant-platform.com">http://oreilly.quant-platform.com</a>

**Derivatives Analytics with Python** |
<a href="http://www.derivatives-analytics-with-python.com" target="_blank">Derivatives Analytics @ Wiley Finance</a>

**Python for Finance** |
<a href="http://shop.oreilly.com/product/0636920032441.do" target="_blank">Python for Finance @ O'Reilly</a>

9560126286dbc6152981abbd08e19fdbde030fd7
ipynb
Jupyter Notebook
ipython/09_Math_Tools.ipynb
portfolioscout/py4fi
9a65df340189ed52037456da221bf66fe89e787f
CNRI-Python
2018-07-10T09:18:23.000Z
2021-12-30T06:35:09.000Z
ipython/09_Math_Tools.ipynb
ioancw/py4fi
bbf7b41d375e4f7b0344bc9b1e97d7910ad1e6ec
CNRI-Python
ipython/09_Math_Tools.ipynb
ioancw/py4fi
bbf7b41d375e4f7b0344bc9b1e97d7910ad1e6ec
CNRI-Python
2018-01-08T01:10:22.000Z
2021-05-26T17:35:35.000Z
Qwen/Qwen-72B
1. YES
2. YES