Numerical Methods in finance¶
Lecture 1: Introduction to Python¶
by: Tomasz Romańczukiewicz
web: th.if.uj.edu.pl/~trom/
rm B-2-03
Outline¶
- Intro
- Python among other programming languages
- Why is the allegedly fast Python so slow?
- Why is it fast?
- IDEs
- First codes
Intro¶
Python among other programming languages¶
Programing languages:
- 700 (wiki) - 8945 (Historical Encyclopaedia of Programming Languages)
- compiled (C/C++, Java, Fortran, Assembler, Pascal ... ) usually faster
program is compiled into machine code and executed afterwards. - interpreted (Bash, Perl, Python, Julia, Octave, Matlab ...) usually much slower
code (strings) are converted (many operations) and executed at the running time
Some of the languages can be "compiled" to speed up (Julia, Python), but usually it requires additional tools and may be complicated
Why is the allegedly fast Python so slow?¶
- Python productivity: usually little code is needed to implement some algorythm so it is fast to write
- Easy to read, the code is usually clean with low $WTF/s$
- It is advertised as fast
- but codes rewritten from standard books run 50-200 times slower than those in C
- Loops (and recursive functions) are interpreted as many times as the loops are executed
- Calculations are encoded by strings
- String operations are much slower and more complicated than calculations on numbers
OK, since it's slow, why is it fast?¶
- It is usually faster than most interpreted languages such as Matlab, Octave
but is slower than some others like Julia (advertised as fast as C as productive as Python) - With additional work some (but not all) problems can be vectorized (for example with numpy) and rewtitten without loops
- Vectorized expresion is interpreted once and all the calculations are done by compiled methods (usually written in highly optimized C) and runs very fast
- Iterations when some value is calculated from previos value, such as $F_n=F_{n-1}+F_{n-2}$ cannot be vectorized and still run 100x slower than C
- Some (partial) solutions:
- Cython - almost Python code translated to C, compiled and run
- numba - functions compiled at run time
It is usually listed as the second best tool for a given problem:
- www: JavaScript is on top but Python is close behind, the list is changing fast
- interpretted language for numerical calculations:
Matlab/Octave beats Python in productivity
Python beats Matlab in speed but Julia is faster - Plotting (my personal list)
- Gnuplot
- Python
- Gnuplot
Example of vectorization¶
Task: Calculate a sume $s=\scriptsize\displaystyle\sum_{n=1}^Nn$ for some large $N$
In [1]:
%%time
N = int(1e7)
s = 0
for n in range(1,N+1):
s += n
In [3]:
%%time
N = int(1e7)
s = np.sum(np.arange(1,N+1))
Python can be¶
- run in iterative shell (just type python3 or ipython3 in a terminal)
and executed on a fly
values are displayed right after execution - edited in some editor (even as simple as vi, vim, gedit, kate etc) and then
run as scripts from commandline or
via some embedded tools - editted dedicated IDEs (Spyder, Pycharm, Pydev etc)
run from within IDE with additional features such as list
output with plots
variable lists with current values - edited, run and presented in notebooks such as this one (Jupyter)
First programs¶
- Beginers should go through these codes as homework: run, change and play with the codes until everything is clear
- write some code by your own
- You can/should study some online tutorials
- Official tutorial
- W3School
- other tutorials and books
- Let me know which worked best for you
- print out some cheat codes
Python as a calculator¶
open iterative console: type python3 or ipython3 in a terminal)
In [4]:
1+2
Out[4]:
In [5]:
0.1+0.2
Out[5]:
In [6]:
0.1+0.2-0.3
Out[6]:
In [7]:
type(1+2)
Out[7]:
In [8]:
type(0.1+0.2)
Out[8]:
In [9]:
1/2
Out[9]:
In [10]:
type(1/2)
Out[10]:
In [11]:
1//2
Out[11]:
In [12]:
type(1//2)
Out[12]:
In [13]:
2**5
Out[13]:
In [14]:
6%4
Out[14]:
In [15]:
help(print)
Iterative console is fun (automatic output) but for serious job we need a proper editor (IDE)
- My choice: Spyder
- many features:
- variable explorer,
- integrated iconsole with possibility to render plots and LaTeX
- nice help
- many features:
- Pycharm
- Inteliji fork for Python
- free but for cool features you have to pay
- Pydev
- but you can use any editor and use its tools to run scripts or run in external terminal
python3 scriptname.py
In [16]:
# Ordinary comment
"""
Multiline documentation
"""
print("Hello world") # print out a string
a = 5 # assign a variable
b = 6
print(a, b) # print out numbers
a, b = 3, 2 # multiple assignment, a tupple
print("a = ", a, " b = " , b)
a, b = b, a # simple swap
print("a = ", a, " b/a = ", b/a, "b//a = ", b//a) # unlike in C,Java etc b/a is not integer,
# but b//a is an integer
print("Formatted output: {:10d}, {:.4f}".format(a, b/a))
print("Conversion from floats to ints:", int(1000*a/b))
print("modulo division:", 15%6)
In [17]:
""" Logic operators """
to_be = True
print("1.", to_be or not to_be)
to_be = False
print("2.", to_be or not to_be)
print("3.", to_be and not to_be)
to_be = 3
print("4.", not to_be)
print("5.", 7 | 2, 7&2, ~7)
print("6.", 1<<3, 3<<3, 25>>2)
print("7.", (7%3!=0), (7%3!=0)and(7%2==0), (666<=777))
In [1]:
import math as mth # including the full math module with prefix mth
print("sin(1) = ", mth.sin(1))
print("cos(1) = ", mth.cos(1))
import numpy as np # including the full numpy module with prefix np, in order not to confuse with math
print("sin(1) = ", np.sin(1))
print("cos(1) = ", np.cos(1))
from math import exp, cosh # including only two functions but without prefix
print("exp(1) = ", exp(1))
print("cosh(1) =", cosh(1))
Exercise:¶
Calculate the kinetic energy using relativistic and newtonian definitions
$$\scriptsize E_{rel} = \frac{mc^2}{\sqrt{1-\frac{v^2}{c^2}}}-mc^2$$$$\scriptsize E_{Newt} = \frac{1}{2}mv^2$$for $v=1 \frac{\text{m}}{\text{s}}$, $c=3\cdot 10^9 \frac{\text{m}}{\text{s}}$, $m=1$kg
In [20]:
from math import sqrt
v, c, m = 1, 3e9, 1
Erel = m*c**2/sqrt(1-v**2/c**2)-m*c**2
ENewt = 0.5*m*v**2
print("E_rel = ", Erel, "\nE_Newt = ", ENewt)
What's wrong?
In [21]:
v**2/c**2
Out[21]:
In [22]:
1-v**2/c**2
Out[22]:
Standard floats (64 bits) can handle only up to 16 digits So the accurate equation reduces to
$$\scriptsize E_{rel} = \frac{mc^2}{\sqrt{1-\frac{v^2}{c^2}}}-mc^2\to\frac{mc^2}{1}-mc^2=0$$which is painfully wrong.
Due to the cancellation errors sometimes it is better to use approximation rather than the full formula
But this approximation is wrong for large values of $v$.
Is there any way to have one formula which works in both cases?
With a bit of magic (short multiplication formulas)
$$\scriptsize E_{rel} = mc^2\left(\frac{c}{\sqrt{c^2-v^2}}-1\right) = mc^2 \frac{\frac{c^2}{{c^2-v^2}}-1}{\frac{c}{\sqrt{c^2-v^2}}+1} = \frac{mv^2}{\sqrt{1-\frac{v^2}{c^2}} + 1-\frac{v^2}{c^2}} $$In [23]:
g = 1-v**2/c**2
E3 = m*v**2/(sqrt(g)+g)
print(E3)
The above formula is mathematically equivalent with the definition but cancellation does not produce errors for small values of $v$.
In [22]:
""" Simple calculations: but what can go wrong? """
from numpy import cosh
print("1. ", 0.3+0.2-0.5) # this is 0
print("2. ", 0.1+0.2-0.3) # this is not 0, beware of round of errors
x = 20
print("3. ", (1-cosh(x))/(1+cosh(x)))
x = 800
print("4. ", (1-cosh(x))/(1+cosh(x))) # this produces overflow error
print("5. ", (1/cosh(x)-1)/(1/cosh(x)+1)) # this produces overflow error
In [27]:
""" List examples """
v = [3, 4, 5, 8, -1] # example array
print ("1:", v)
print ("2:", v[0], v[4], v[-2]) # elements: [(3), (4), 5, (8), -1]
# indexes: 0 1 2 3 4
# negative ind: -5 -4 -3 -2 -1
v.append(7) # adding a new element
print("3:", v)
v.insert(2, 6) # inserting at certain position
print("4:", v)
print("5:", v[2:4]) # slices [3, 4, (6, 5), 8, -1, 7]
v.sort() # hmm, what could that be?
print("6:", v)
In [43]:
w = list(range(1, 15, 3)) # each 3rd integer number from [1,15)
print (" 7:", w)
print(" 8:", w[2:5]+v[0:4:2]) # join list slices [1, 4, (7, 10, 13)] [(-1), 3, (4), 5, 6, 7, 8]
# indices: 0 1 2 3 4 0 1 2 3 4 5 6
print(" 9:", 2*w) # the list is repeated twice
w.reverse()
print("10:", w)
print("11:", w.index(7))
w.append("The last in line")
print(w)
print("12:", w.pop())
print("13:", w)
In [26]:
x = (1, 3, 6.62354, "Hello", "world") # example of a tuple
print(x)
print(x[3], x[4])
#x[1] = 2 # an error occurs, tuple as not mutable
y = {2,3,4,2,3, "hello", "hello"} # example of a set
print(y)
# print(y[2]) # error: 'set' object does not support indexing
z = {'a': 5, 'b': 3.1415, 'c': "some text"} # example of a dictionary
print(z)
print(z['c'])
In [27]:
""" Flow control """
from math import sqrt
a = 1; b = 2; c = -1;
Δ = b**2-4*a*c
print ("Δ = ", Δ)
if Δ <0:
print ("Δ is negative, no real solutions")
elif Δ==0: # note: '==' means comparison, '=' is assignment
print ("Δ is equal to 0, one real solution x =", -b/(2*a))
else:
print ("Δ is positive, two real solutions")
x1 = (-b+sqrt(Δ))/(2*a)
x2 = (-b-sqrt(Δ))/(2*a)
print ("x1 = ", x1, "x2 =", x2)
Excercise:¶
- Run the code to see all the possibilities.
- Check $a=c=10^{-4},\;b=2\cdot10^4$. Solve it analytically as well.
- Extend the program to include the case for $a=0$ (linear equation) and consider all posible cases ($b=0, c\neq 0$ and $b=0, c=0$)
In [28]:
""" Loops """
for n in range(2, 15, 3):
print (n)
In [29]:
F = [1, 1]
for n in range(2, 10):
F.append(F[n-2]+F[n-1])
print ("Fibbonaci series: ", F)
In [30]:
a = 120
divisors =[]
for k in range(2,a+1):
if a%k==0:
divisors.append(k)
print ("Divisors of ", a, ": ", divisors)
In [46]:
""" control and some more """
for n in range(10):
if (n%2==0): continue # skip the rest and go to the beginning of the loop
if (n==7): break # finish the loop when n is equal to 7
print(n)
In [32]:
lst=["This", 'is', "a", 3.1415, "complicated list"]
for idx, element in enumerate(lst):
print("index=", idx, "value=", element)
In [33]:
squares = [k**2 for k in range(10)]
print (squares)
In [56]:
"""Iterations"""
x, k = 0.3, 4
for n in range(20):
x = k*x*(1-x)
print(x)
Consider a logistic map $\scriptsize x\ni[0,1]\to f(x)=kx(1-x)\in[0,1]$ for $\scriptsize k\in[0,4]$.
The iteration
$$\scriptsize x_{n+1}=kx_n(1-x_n)$$has two fixed points (solution to $\scriptsize x=f(x)$): $\scriptsize x=0$ and $\scriptsize x=1-\frac1k$.
Q: What happens when $\scriptsize k=0.5$, $\scriptsize k=2.5$, $\scriptsize k=3.1$, $\scriptsize k=3.5$, $\scriptsize k=4$? And why?¶
Around a fixed point $\scriptsize x_*=f(x_*)$:
$$\scriptsize x_n=x_*+h_n$$$$\scriptsize h_{n+1}=x_{n+1}-x_*=f(x_*+h_n)-x_*$$Linear stability
$$\scriptsize h_{n+1}\approx f(x_*)+h_nf'(x_*)-x_*+\mathcal{O}(h_n^2)$$Neglecting nonlinear terms we have a geometric series:
$$\scriptsize h_{n+1}=h_nf'(x_*)$$Fixed point $x_*$ of an iteration $x_{n+1}=f(x_n)$ is linearly stable when $|f'(x_*)|<1$ and is unstable when $|f'(x_*)|>1$. Moreover the smaller $|f'(x_*)|$ the faster the convergence.
$$\scriptsize f'(0)=k\qquad f'\left(1-\frac1k\right)=2-k$$- $k<1$ $x=0$ is the only stable point,
- for $1<k<2$ $x=1-\frac1k$ is the only stable piont,
- for $k>3$ there are no stable points
- but there are stable 2-cycles, 4-cycles etc $x_{n+2}=x_n$
- for certain values of $k$ there is pure chaos
It is interesting from the point of view of a dynamical systems but very bad as a method for solving equations. But we can influence the stability
$$\scriptsize \biggl.x = f(x)\;\biggr|\;+\omega x$$$$\scriptsize (1+\omega)x = f(x)+\omega x$$$$\scriptsize x = \frac{f(x)+\omega x}{1+\omega}\equiv\tilde f(x)$$Parameter $\omega$ can control the stability ad convergence rate
$$\scriptsize \tilde f'(0)=\frac{k+\omega}{1+\omega}\qquad \tilde f'\left(1-\frac1k\right)=\frac{2-k+\omega}{1+\omega}$$In [35]:
"""Increase stabililty of the iteration method """
x, k, ω = 0.5, 4, 1.5
for n in range(20):
x = ( k*x*(1-x) + ω*x) / (1+ω)
print(x)
Another simple iteration:¶
$$\scriptsize I_n=\frac{1}{e}\int_0^1e^xx^n\;dx$$generates the following conditions (decreasing but positive sequence)
$$\scriptsize 0<I_n<I_{n-1}$$and
$$\scriptsize I_0=1-\frac{1}{e}\qquad I_n=1-nI_{n-1}$$In [53]:
import numpy as np
I = 1-1/np.e
for n in range(1,25):
if n>15: print("{:2d}\t{:15.8f}".format(n, I))
I = 1-n*I
In [37]:
""" Functions """
def Function():
print("This is a simple function. Just some code")
Function(); Function() # function calls
In [38]:
def Sum(a, b): # arguments
return a*b; # returned value
result = Sum(3, 222)
print(result, Sum(3,4))
In [39]:
def g(a, b):
return (a+b, a*b, a/b, a**b) # return a tuple of four values
s, m, d, p = g(3, 4)
print(s, m, d, p)
tup = g(3, 4)
print(tup)
Implementation of the recursive definition $$\scriptsize n!=\left\{ \begin{array}{cl} 1&\text{ for } n=0\\ n\,(n-1)!&\text{ otherwise } \end{array} \right.$$
In [44]:
""" Recursive functions """
def Factorial(n):
if n<1: return 1
else: return n*Factorial(n-1)
Factorial(4)
Out[44]:
In [45]:
""" Recursive functions """
def Factorial(n):
print ("starting the calculation of ", n, "!")
if n<1: result = 1
else: result = n*Factorial(n-1)
print ("finished the calculation of ", n, "! = ", result)
return result
Factorial(4);
Fibonacci series: $F_1=F_2=1$ and $F_{n}=F_{n-1}+F_{n-2}$ for $n>2$
In [6]:
def Fib(n):
if n>2: return Fib(n-1)+Fib(n-2)
else: return 1
Fib(10)
Out[6]:
In [62]:
"""Very bad Fibbonacci"""
counter = 0
def Fib(n):
global counter # global counter
counter+=1
if n<3: return 1
else: return Fib(n-1)+Fib(n-2)
print("Fib(10) = ", Fib(10))
print("counter = ", counter)
In [70]:
""" Functions as arguments """
import math as mt
def make_table(f):
print("==========",f.__name__,"============")
for n in range(10):
#print(f"{n*0.1:.2f}\t\t{f(n*0.1): .8f}")
print("{:.2f}\t\t{: .8f}". format(n*0.1, f(n*0.1)))
make_table(mt.sin) # another function as an argument
In [71]:
make_table(lambda x: x**2-1) # anonymous lambda function
In [74]:
"""Functions with default and key arguments """
def f1(a, b=3, c=4):
print("f1:", a, b, c)
f1(3, 4)
f1(2, c=7)
In [75]:
def f2(a, *b): # additional arguments are used as a tuple
print("f2:", a, b)
f2("Hello")
f2("Hello", "world", 777)
In [73]:
def f3(**k): # all arguments are a dictionary
print("f3:", k)
for n in k:
print("f3: ", n," = ", k[n])
f3(a=1, b=2, word="Hello")
Summary of Python¶
- Python is a very intuitive and powerfull language.
- Basit types are strings, integers and floating-point numbers.
- They can be encapsulated in lists [], tuples () and sets {}.
- Lists are mutable, tuples and sets are not.
- Division / of two numbers always give a floating point number.
- Integer division: //, power: **.
- if ..else statement controls the flow of the programs (decision making).
- for loop can loop thourgh the lists, butis't slow.
- Functions can be defined to use fragments of codes efficiently.
- Anonymous lambda functions can be defined, when small function can be used as an argument.
Summary of the numerical part¶
- All floating point calculations generate errors.
- Even simple equations can give huge relative errors (cancelation problem).
- Sometimes formulea can be rewritten in more numerical friendly way.
- Sometimes approximations are better than exact formulae.
- Badly written program can lead to extreme complexity (recursive Fibbonacci)
- Iterative methods can have lots of problems:
- stability issues,
- two-,four- or more cycles,
- Sometimes modification can change convergence and stability.
- Some error can accumulate very quickly.
Calculations in Python/numpy¶
In [76]:
import numpy as np # np is used as a namespace
x_list = list(range(5))
print("list: ", x_list)
x = np.arange(5)
print("np.array:", x)
print("2*x_list:", 2*x_list) # lists are joined
x = np.arange(5)
print("2*x: ", 2*x) # arrays are calculated elementwise
print("cubes: ", x**3)
In [81]:
x = np.linspace(0, np.pi, 8) # horizontal vector
print ("shape of x:", np.shape(x))
print ("\ntable of sins: ", np.sin(x))
y = np.array( [x, np.sin(x), np.cos(x), np.sin(x)**2+np.cos(x)**2 ] )
# 4 rows of horizontal vectors
print ("\nshape of y", np.shape(y))
print (y.transpose()) # transpose for nice output
In [80]:
sums=np.zeros(4)
print ("\nsums after initiation:", sums)
for k in range(4):
sums[k] = np.sum(y[:,k])
print ("Calculated sums:", sums)
Recall the example: $\scriptsize\displaystyle S = \sum_{n=1}^{N}n$.
Numpy can be much much faster than standard Python loops
In [1]:
%%time
N=int(1e7)
s = 0
for n in range(1,N+1):
s += n
In [7]:
%%time
N=int(1e7)
s = np.sum(np.arange(1,N+1))
In [8]:
"""Broadcasting examples"""
a = np.array([1, 2, 3])
3*a # 3*a_i
Out[8]:
In [9]:
a*a # a_i^2
Out[9]:
In [18]:
b = np.ones([2,3])
2*b
Out[18]:
In [19]:
a*b
Out[19]:
In [20]:
b*a
Out[20]:
In [106]:
import matplotlib.pyplot as plt # ploting library
def make_plot():
x = np.linspace(0, 2*np.pi, 1000)
y = np.array( [ np.sin(x), np.cos(x), np.exp(-0.5*x)*np.sin(5*x) ] ).transpose()
plt.figure(figsize=(8,4), dpi=120)
plt.plot(x, y)
plt.legend(['$\sin(x)$', '$\cos(x)$', '$e^{-x/2}\sin(5x)$'])
plt.title('Example plots')
plt.xlabel('$x$')
plt.ylabel("$y$")
plt.grid(True)
plt.show()
In [107]:
make_plot()
In [85]:
"""Sympy example"""
import sympy as sp
x, k = sp.var('x k')
f = k*x*(1-x)
print("Expression:", f)
eq1 = sp.Eq(f-x)
sols = sp.solve(eq1,x)
print("Fixed points:", sols)
In [86]:
df = sp.diff(f,x)
print("derivative:", df)
for (n,s) in enumerate(sols):
print("f'(x_",n,") = ", sp.simplify(df.subs(x, s)))
In [101]:
g=f.subs(x,f) # g(x) = f(f(x))
print(g, "\n")
two_cycle = sp.solve(sp.Eq(g-x),x)
print(*two_cycle, "\n", sep="\n")
print( sp.latex(two_cycle[2]))
print(sp.latex(two_cycle[3]))
In [105]: