<br>

# Python for Finance (2nd ed.)

**Mastering Data-Driven Finance**

&copy; Dr. Yves J. Hilpisch | The Python Quants GmbH



# Mathematical Tools

## Approximation


```python
import numpy as np
from pylab import plt, mpl
```


```python
plt.style.use('seaborn')
mpl.rcParams['font.family'] = 'serif'
%config InlineBackend.figure_format = 'svg'
```


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


```python
def create_plot(x, y, styles, labels, axlabels):
    plt.figure(figsize=(10, 6))
    for i in range(len(x)):
        plt.plot(x[i], y[i], styles[i], label=labels[i])
        plt.xlabel(axlabels[0])
        plt.ylabel(axlabels[1])
    plt.legend(loc=0)
```


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


```python
create_plot([x], [f(x)], ['b'], ['f(x)'], ['x', 'f(x)'])
# plt.savefig('../../images/ch11/math_01.png');
```


    

    


### Regression

#### Monomials as Basis Functions


```python
res = np.polyfit(x, f(x), deg=1, full=True)  
```


```python
res  
```




    (array([ 4.28841952e-01, -4.31775426e-17]),
     array([21.03238686]),
     2,
     array([1., 1.]),
     1.1102230246251565e-14)




```python
ry = np.polyval(res[0], x)  
```


```python
create_plot([x, x], [f(x), ry], ['b', 'r.'],
            ['f(x)', 'regression'], ['x', 'f(x)'])
# plt.savefig('../../images/ch11/math_02.png');
```


    

    



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


```python
create_plot([x, x], [f(x), ry], ['b', 'r.'],
            ['f(x)', 'regression'], ['x', 'f(x)'])
# plt.savefig('../../images/ch11/math_03.png');
```


    

    



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


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




    False




```python
np.mean((f(x) - ry) ** 2)  
```




    0.0017769134759517873




```python
create_plot([x, x], [f(x), ry], ['b', 'r.'],
            ['f(x)', 'regression'], ['x', 'f(x)'])
# plt.savefig('../../images/ch11/math_04.png');
```


    

    


#### 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), rcond=None)[0]  
```


```python
reg.round(4)  
```




    array([ 0.    ,  0.5628, -0.    , -0.0054])




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


```python
create_plot([x, x], [f(x), ry], ['b', 'r.'],
            ['f(x)', 'regression'], ['x', 'f(x)'])
# plt.savefig('../../images/ch11/math_05.png');
```


    

    



```python
matrix[3, :] = np.sin(x)  
```


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


```python
reg.round(4)  
```




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




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


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




    True




```python
np.mean((f(x) - ry) ** 2)  
```




    2.099263639382087e-31




```python
create_plot([x, x], [f(x), ry], ['b', 'r.'],
            ['f(x)', 'regression'], ['x', 'f(x)'])
# plt.savefig('../../images/ch11/math_06.png');
```


    

    


#### 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
create_plot([x, x], [f(x), ry], ['b', 'r.'],
            ['f(x)', 'regression'], ['x', 'f(x)'])
# plt.savefig('../../images/ch11/math_07.png');
```


    

    


#### 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))  
```

    [ 0.2   3.85  4.64  5.67 -2.    3.41  2.8  -4.69  3.21 -5.26]
    [ 0.3   1.27  1.32  2.26 -1.91  1.44  1.73 -1.35  1.54 -1.77]



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


```python
create_plot([xu, xu], [yu, ry], ['b.', 'ro'],
            ['f(x)', 'regression'], ['x', 'f(x)'])
# plt.savefig('../../images/ch11/math_08.png');
```


    

    


#### Multiple Dimensions


```python
def fm(p):
    x, y = p
    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)  
```


```python
Z = fm((X, Y))
x = X.flatten()  
y = Y.flatten()  
```


```python
from mpl_toolkits.mplot3d import Axes3D  
```


```python
fig = plt.figure(figsize=(10, 6))
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, Z, rstride=2, cstride=2,
                       cmap='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)
# plt.savefig('../../images/ch11/math_09.png');
```

    <ipython-input-40-4729664be1fa>:2: MatplotlibDeprecationWarning: Calling gca() with keyword arguments was deprecated in Matplotlib 3.4. Starting two minor releases later, gca() will take no keyword arguments. The gca() function should only be used to get the current axes, or if no axes exist, create new axes with default keyword arguments. To create a new axes with non-default arguments, use plt.axes() or plt.subplot().
      ax = fig.gca(projection='3d')





    <matplotlib.colorbar.Colorbar at 0x113e452b0>




    

    



```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
reg = np.linalg.lstsq(matrix, fm((x, y)), rcond=None)[0]
```


```python
RZ = np.dot(matrix, reg).reshape((20, 20))  
```


```python
fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, figsize=(10, 6))
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, ax=ax)
# plt.savefig('../../images/ch11/math_10.png');
```




    <matplotlib.colorbar.Colorbar at 0x113fbf730>




    

    


### 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
np.allclose(f(x), iy)  
```




    True




```python
create_plot([x, x], [f(x), iy], ['b', 'ro'],
            ['f(x)', 'interpolation'], ['x', 'f(x)'])
# plt.savefig('../../images/ch11/math_11.png');
```


    

    



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


```python
create_plot([xd, xd], [f(xd), iyd], ['b', 'ro'],
            ['f(x)', 'interpolation'], ['x', 'f(x)'])
# plt.savefig('../../images/ch11/math_12.png');
```


    

    



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


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




    False




```python
np.mean((f(xd) - iyd) ** 2)  
```




    1.1349319851457788e-08




```python
create_plot([xd, xd], [f(xd), iyd], ['b', 'ro'],
            ['f(x)', 'interpolation'], ['x', 'f(x)'])
# plt.savefig('../../images/ch11/math_13.png');
```


    

    


## Convex Optimization


```python
def fm(p):
    x, y = p
    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, ax = plt.subplots(subplot_kw={"projection": "3d"}, figsize=(10, 6))
surf = ax.plot_surface(X, Y, Z, rstride=2, cstride=2,
                       cmap='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)
# plt.savefig('../../images/ch11/math_14.png');
```




    <matplotlib.colorbar.Colorbar at 0x11520bac0>




    

    


### Global Optimization


```python
import scipy.optimize as sco  
```


```python
def fo(p):
    x, y = p
    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
sco.brute(fo, ((-10, 10.1, 5), (-10, 10.1, 5)), finish=None)  
```

    -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 = sco.brute(fo, ((-10, 10.1, 0.1), (-10, 10.1, 0.1)), finish=None)
```


```python
opt1
```




    array([-1.4, -1.4])




```python
fm(opt1)
```




    -1.7748994599769203



### Local Optimization


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

     -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.



```python
opt2
```




    array([-1.42702972, -1.42876755])




```python
fm(opt2)
```




    -1.7757246992239009




```python
output = False
sco.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
import math
```


```python
def Eu(p):  
    s, b = p
    return -(0.5 * math.sqrt(s * 15 + b * 5) +
             0.5 * math.sqrt(s * 5 + b * 12))
```


```python
cons = ({'type': 'ineq',
         'fun': lambda p: 100 - p[0] * 10 - p[1] * 10})  
```


```python
bnds = ((0, 1000), (0, 1000))  
```


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


```python
result
```




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




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




    array([8.02547122, 1.97452878])




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




    9.700883611487832




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




    100.0



## Integration


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


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


```python
x = np.linspace(0, 10)
y = f(x)
a = 0.5  
b = 9.5  
Ix = np.linspace(a, b)  
Iy = f(Ix)  #  <4>
```


```python
from matplotlib.patches import Polygon
```


```python
fig, ax = plt.subplots(figsize=(10, 6))
plt.plot(x, y, 'b', linewidth=2)
plt.ylim(bottom=0)
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)
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)]);
# plt.savefig('../../images/ch11/math_15.png');
```


    

    


### Numerical Integration


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




    24.3669959670846




```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.37496418455075



### 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.mean(f(x)) * (b - a))
```

    24.804762279331463
    26.522918898332378
    26.265547519223976
    26.02770339943824
    24.99954181440844
    23.881810141621663
    23.527912274843253
    23.507857658961207
    23.67236746066989
    23.679410416062886
    24.424401707879305
    24.239005346819056
    24.115396924962802
    24.424191987566726
    23.924933080533783
    24.19484212027875
    24.117348378249833
    24.100690929662274
    23.76905109847816


## Symbolic Computation

### Basics


```python
import sympy as sy
```


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


```python
type(x)
```




    sympy.core.symbol.Symbol




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




$\displaystyle \sqrt{x}$




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




$\displaystyle \sqrt{x} - 13$




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


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




$\displaystyle 1.5 x^{2} + 4.5$




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


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

         2      
    1.5*x  + 4.5



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

      ___      
    \/ x  + 0.5



```python
%time pi_str = str(sy.N(sy.pi, 400000))  
```

    CPU times: user 170 ms, sys: 10.4 ms, total: 181 ms
    Wall time: 217 ms



```python
pi_str[:42]  
```




    '3.1415926535897932384626433832795028841971'




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




    '8245672736856312185020980470362464176199'




```python
%time pi_str.find('061072')  
```

    CPU times: user 81 µs, sys: 0 ns, total: 81 µs
    Wall time: 83 µs





    80847



### 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 and Differentiation


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


```python
I = sy.Integral(sy.sin(x) + 0.5 * x, (x, a, b))  
```


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

      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  
```




    24.3747547180867




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


```python
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.3747547180867



### 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  
```




    -1.42755177876459




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




    -1.42755177876459




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




    -1.77572565314742




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




    1.74632928225285




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




    1.74632928225285




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




    2.27423381055640



<br>

<a href="http://tpq.io" target="_blank">http://tpq.io</a> | <a href="http://twitter.com/dyjh" target="_blank">@dyjh</a> | <a href="mailto:training@tpq.io">training@tpq.io</a>

b4cf74673627e7d2fb2622afac7a7fab029bb8c9
ipynb
Jupyter Notebook
code/ch11/11_math_tools.ipynb
Lyricsy/py4fi2nd
4dba185b193a618aab618305d45d62218f2bab43
CNRI-Python
2018-12-01T12:31:56.000Z
2022-03-31T03:42:26.000Z
code/ch11/11_math_tools.ipynb
claire-caicai/py4fi2nd
ae94e80d45e277eb51af21ce6c1765c07c5dbcbe
CNRI-Python
2018-12-03T07:46:25.000Z
2021-12-17T14:22:46.000Z
code/ch11/11_math_tools.ipynb
claire-caicai/py4fi2nd
ae94e80d45e277eb51af21ce6c1765c07c5dbcbe
CNRI-Python
2018-12-03T01:46:56.000Z
2022-03-27T09:43:49.000Z
Qwen/Qwen-72B
1. YES
2. YES