<jupyter_start><jupyter_text>SciPy Author: <NAME>SciPy is a collection of mathematical algorithms and convenience functions. In this this notebook there are just a few examples of the features that are most important to us. But if you want to see all that SciPy has to offer, have a look at the [official documentation](https://docs.scipy.org/doc/scipy/reference/).Since SciPy has several subpackages, it is commom practice to import just the one we are going to use, as you'll in the following examples.<jupyter_code>import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Definite IntegralsThe function `quad` is provided to integrate a function of one variable between two points. This functions has 2 outputs, the first one is the computed integral value and the second is an estimate of the absolute error.<jupyter_code>import scipy.integrate as integrate

def my_func(x):
    return x**2

integrate.quad(my_func, 0, 2)<jupyter_output><empty_output><jupyter_text>The `quad` functions also allows for infinite limits.$$\int_{-\infty}^{\infty} e^{-x^{2}}dx$$<jupyter_code>def my_func(x):
    return np.exp(-x**2)

integrate.quad(my_func, -np.inf, np.inf)<jupyter_output><empty_output><jupyter_text>SciPy's `integrate` library also has functions for double and triple integrals. Check them out in the documentations.  OptimizationThe `scipy.optimize` package provides several commonly used optimization algorithms. Here we are going to use just one to illustrate.Consider that you have 3 assets available. Their expected returns, risks (standard-deviations) and betas are on the table bellow and $\rho$ is the correlation matrix of the returns.| Asset | Return | Risk | Beta ||-------|--------|------|------||A      |3%      | 10%  | 0.5  ||B      |3.5%    | 11%  | 1.2  ||C      |5%      | 15%  | 1.8  |$$\rho = \begin{bmatrix}1 & 0.3 & -0.6 \\0.3 & 1 & 0 \\-0.6 & 0 & 1 \end{bmatrix}$$Use the `minimize` function to find the weights of each asset that maximizes it's Sharpe index.<jupyter_code>retu = np.array([0.03, 0.035, 0.05])
risk = np.array([0.10, 0.11, 0.15])
beta = np.array([0.5, 1.2, 1.8])

corr = np.array([[1, 0.3, -0.6], 
                 [0.3, 1, 0],
                 [-0.6, 0, 1]])

def port_return(w):
    return retu.dot(w)

def port_risk(w):
    covar = np.diag(risk).dot(corr).dot(np.diag(risk))
    return (w.dot(covar).dot(w))**0.5

def port_sharpe(w):
    return -1*(port_return(w) / port_risk(w))   # The -1 is because we want to MINIMIZE the negative of the Sharpe

def port_weight(w):
    return w.sum()<jupyter_output><empty_output><jupyter_text>When declaring an optimization problem with inequality restrictions, they have the form of:$$\begin{align*}\min_{w} & f\left(w\right)\\s.t. & g\left(w\right)\geq0\end{align*}$$<jupyter_code>from scipy.optimize import minimize

eq_cons = {'type': 'eq',
           'fun' : lambda w: port_weight(w) - 1}

w0 = np.array([1, 0, 0])

res = minimize(port_sharpe, w0, method='SLSQP', constraints=eq_cons, options={'ftol': 1e-9, 'disp': True})
res.x
res.x.sum()
-1*res.fun<jupyter_output><empty_output><jupyter_text>InterpolationThere are several general interpolation facilities available in SciPy, for data in 1, 2, and higher dimensions. The `interp1d` funtions grabs data points and returns a function. The default interpolation method is the linear interpolation, but there are several to choose from.<jupyter_code>from scipy.interpolate import interp1d

x = np.linspace(0, 10, num=11, endpoint=True)
y = np.cos(-x**2/9.0)
f1 = interp1d(x, y)  # linear is the default
f2 = interp1d(x, y, kind='cubic')  # cubic splines
f3 = interp1d(x, y, kind='nearest')  # grab the nearest value
f4 = interp1d(x, y, kind='previous')  # hold last value
f5 = interp1d(x, y, kind='next')  # grab the next value<jupyter_output><empty_output><jupyter_text>Now that we have the interpolated function, lets generate more points for the x axis and plot the different methods<jupyter_code>xnew = np.linspace(0, 10, num=41, endpoint=True)
xnew
plt.plot(x, y, 'o', xnew, f1(xnew), '-', xnew, f2(xnew), '--')
plt.legend(['data', 'linear', 'cubic'], loc='best')
plt.show()
plt.plot(x, y, 'o', xnew, f3(xnew), '-', xnew, f4(xnew), '--', xnew, f5(xnew), ':')
plt.legend(['data', 'nearest', 'previous', 'next'], loc='best')
plt.show()<jupyter_output><empty_output><jupyter_text>The `interpolate` sublibrary also has interpolation methods for multivariate data and has **integration with pandas**. Have a look at the documentation.  Linear Algebra (again)`scipy.linalg` contains all the functions in `numpy.linalg` plus some more advanced ones.<jupyter_code>from scipy import linalg as la

A = np.array([[1,3,5],[2,5,1],[2,3,8]])
la.inv(A)<jupyter_output><empty_output><jupyter_text>Matrix and vector **norms** can also be computed with SciPy. A wide range of norm definitions are available using different parameters to the order argument of `linalg.norm`.<jupyter_code>A = np.array([[1, 2], [3, 4]])
print(la.norm(A))  # frobenius norm is the default.
print(la.norm(A, 1)) # L1 norm (max column sum)
print(la.norm(A, np.inf)) # L inf norm (max row sum)<jupyter_output>5.477225575051661
6.0
7.0<jupyter_text>Some more advanced matrix decompositions are also available, like the **Schur Decomposition**<jupyter_code>la.schur(A)<jupyter_output><empty_output><jupyter_text>Some notable matrices can also be created, like block **diagonal matrices**.<jupyter_code>A = np.array([[1, 0],
              [0, 1]])

B = np.array([[3, 4, 5],
              [6, 7, 8]])

C = np.array([[7]])

la.block_diag(A, B, C)<jupyter_output><empty_output><jupyter_text>Solving Linear Systems$$\begin{align}x+3y+5 & =10\\2x+5y+z & =8\\2x+3y+8z & =3\end{align}$$The system above can be written with matrix notation as $AX=B$ and we know we can find the solution by doing $X=A^{-1}B$, but inverting a matrix is computationally expensive. When solving big linear system it is advised to use the `solve` method.<jupyter_code>A = np.array([[1, 3, 5], [2, 5, 1], [2, 3, 8]])
B = np.array([[10], [8], [3]])<jupyter_output><empty_output><jupyter_text>Lets check the time that it takes to solve the system in both ways...<jupyter_code>la.inv(A).dot(B)
la.solve(A, B)<jupyter_output><empty_output><jupyter_text>let's try with a bigger matrix<jupyter_code>import numpy.random as rnd
A = rnd.random((1000, 1000))
B = rnd.random((1000, 1))
%%timeit
la.inv(A).dot(B)
%%timeit
la.solve(A, B)<jupyter_output>14.9 ms ± 495 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)<jupyter_text>SciPy Author: <NAME>SciPy is a collection of mathematical algorithms and convenience functions. In this this notebook there are just a few examples of the features that are most important to us. But if you want to see all that SciPy has to offer, have a look at the [official documentation](https://docs.scipy.org/doc/scipy/reference/).Since SciPy has several subpackages, it is commom practice to import just the one we are going to use, as you'll in the following examples.<jupyter_code>import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Definite IntegralsThe function `quad` is provided to integrate a function of one variable between two points. This functions has 2 outputs, the first one is the computed integral value and the second is an estimate of the absolute error.<jupyter_code>import scipy.integrate as integrate

def my_func(x):
    return x**2

integrate.quad(my_func, 0, 2)<jupyter_output><empty_output><jupyter_text>The `quad` functions also allows for infinite limits.$$\int_{-\infty}^{\infty} e^{-x^{2}}dx$$<jupyter_code>def my_func(x):
    return np.exp(-x**2)

integrate.quad(my_func, -np.inf, np.inf)<jupyter_output><empty_output><jupyter_text>SciPy's `integrate` library also has functions for double and triple integrals. Check them out in the documentations.  OptimizationThe `scipy.optimize` package provides several commonly used optimization algorithms. Here we are going to use just one to illustrate.Consider that you have 3 assets available. Their expected returns, risks (standard-deviations) and betas are on the table bellow and $\rho$ is the correlation matrix of the returns.| Asset | Return | Risk | Beta ||-------|--------|------|------||A      |3%      | 10%  | 0.5  ||B      |3.5%    | 11%  | 1.2  ||C      |5%      | 15%  | 1.8  |$$\rho = \begin{bmatrix}1 & 0.3 & -0.6 \\0.3 & 1 & 0 \\-0.6 & 0 & 1 \end{bmatrix}$$Use the `minimize` function to find the weights of each asset that maximizes it's Sharpe index.<jupyter_code>retu = np.array([0.03, 0.035, 0.05])
risk = np.array([0.10, 0.11, 0.15])
beta = np.array([0.5, 1.2, 1.8])

corr = np.array([[1, 0.3, -0.6], 
                 [0.3, 1, 0],
                 [-0.6, 0, 1]])

def port_return(w):
    return retu.dot(w)

def port_risk(w):
    covar = np.diag(risk).dot(corr).dot(np.diag(risk))
    return (w.dot(covar).dot(w))**0.5

def port_sharpe(w):
    return -1*(port_return(w) / port_risk(w))   # The -1 is because we want to MINIMIZE the negative of the Sharpe

def port_weight(w):
    return w.sum()<jupyter_output><empty_output><jupyter_text>When declaring an optimization problem with inequality restrictions, they have the form of:$$\begin{align*}\min_{w} & f\left(w\right)\\s.t. & g\left(w\right)\geq0\end{align*}$$<jupyter_code>from scipy.optimize import minimize

eq_cons = {'type': 'eq',
           'fun' : lambda w: port_weight(w) - 1}

w0 = np.array([1, 0, 0])

res = minimize(port_sharpe, w0, method='SLSQP', constraints=eq_cons, options={'ftol': 1e-9, 'disp': True})
res.x
res.x.sum()
-1*res.fun<jupyter_output><empty_output><jupyter_text>InterpolationThere are several general interpolation facilities available in SciPy, for data in 1, 2, and higher dimensions. The `interp1d` funtions grabs data points and returns a function. The default interpolation method is the linear interpolation, but there are several to choose from.<jupyter_code>from scipy.interpolate import interp1d

x = np.linspace(0, 10, num=11, endpoint=True)
y = np.cos(-x**2/9.0)
f1 = interp1d(x, y)  # linear is the default
f2 = interp1d(x, y, kind='cubic')  # cubic splines
f3 = interp1d(x, y, kind='nearest')  # grab the nearest value
f4 = interp1d(x, y, kind='previous')  # hold last value
f5 = interp1d(x, y, kind='next')  # grab the next value<jupyter_output><empty_output><jupyter_text>Now that we have the interpolated function, lets generate more points for the x axis and plot the different methods<jupyter_code>xnew = np.linspace(0, 10, num=41, endpoint=True)
xnew
plt.plot(x, y, 'o', xnew, f1(xnew), '-', xnew, f2(xnew), '--')
plt.legend(['data', 'linear', 'cubic'], loc='best')
plt.show()
plt.plot(x, y, 'o', xnew, f3(xnew), '-', xnew, f4(xnew), '--', xnew, f5(xnew), ':')
plt.legend(['data', 'nearest', 'previous', 'next'], loc='best')
plt.show()<jupyter_output><empty_output><jupyter_text>The `interpolate` sublibrary also has interpolation methods for multivariate data and has **integration with pandas**. Have a look at the documentation.  Linear Algebra (again)`scipy.linalg` contains all the functions in `numpy.linalg` plus some more advanced ones.<jupyter_code>from scipy import linalg as la

A = np.array([[1,3,5],[2,5,1],[2,3,8]])
la.inv(A)<jupyter_output><empty_output><jupyter_text>Matrix and vector **norms** can also be computed with SciPy. A wide range of norm definitions are available using different parameters to the order argument of `linalg.norm`.<jupyter_code>A = np.array([[1, 2], [3, 4]])
print(la.norm(A))  # frobenius norm is the default.
print(la.norm(A, 1)) # L1 norm (max column sum)
print(la.norm(A, np.inf)) # L inf norm (max row sum)<jupyter_output>5.477225575051661
6.0
7.0<jupyter_text>Some more advanced matrix decompositions are also available, like the **Schur Decomposition**<jupyter_code>la.schur(A)<jupyter_output><empty_output><jupyter_text>Some notable matrices can also be created, like block **diagonal matrices**.<jupyter_code>A = np.array([[1, 0],
              [0, 1]])

B = np.array([[3, 4, 5],
              [6, 7, 8]])

C = np.array([[7]])

la.block_diag(A, B, C)<jupyter_output><empty_output><jupyter_text>Solving Linear Systems$$\begin{align}x+3y+5 & =10\\2x+5y+z & =8\\2x+3y+8z & =3\end{align}$$The system above can be written with matrix notation as $AX=B$ and we know we can find the solution by doing $X=A^{-1}B$, but inverting a matrix is computationally expensive. When solving big linear system it is advised to use the `solve` method.<jupyter_code>A = np.array([[1, 3, 5], [2, 5, 1], [2, 3, 8]])
B = np.array([[10], [8], [3]])<jupyter_output><empty_output><jupyter_text>Lets check the time that it takes to solve the system in both ways...<jupyter_code>la.inv(A).dot(B)
la.solve(A, B)<jupyter_output><empty_output><jupyter_text>let's try with a bigger matrix<jupyter_code>import numpy.random as rnd
A = rnd.random((1000, 1000))
B = rnd.random((1000, 1))
%%timeit
la.inv(A).dot(B)
%%timeit
la.solve(A, B)<jupyter_output>14.9 ms ± 495 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
{'fraction_non_alphanumeric': 0.14571056860794368, 'fraction_numerical': 0.03275426647279406, 'mean_word_length': 4.6498054474708175, 'pattern_counts': {'":': 0, '<': 111, '<?xml version=': 0, '>': 111, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '28', 'dir': 'jupyter-structured-clean-dedup', 'id': '344284', 'license': 'MIT', 'n_tokens_mistral': 4926, 'n_tokens_neox': 4608, 'n_words': 1447, 'path': 'fhnotebooks/Introduction to Python/Section 03 - SciPy.ipynb', 'repo_name': 'fernandotm1/FinanceHub', 'size': '13087'}