Today I will try to benchmark the execution speed of several programming languages on a Monte Carlo example. This benchmark involves VBA, C++, C#, Python, Cython and Numpy vectorization. I will try to add progressively other programming languages so that this article will be more thorough.
Execution environment
All the chunks of code have been executed on the same computer so that the benchmark has sense. It’s a MacBook Pro with 8go RAM, i5 Dual-Core 2.4GHZ, and SSD 256go.
Python code is launch from the console, C# code in Visual Studio 2015, VBA in Excel 2016 and C++ in CodeBlocks with GCC compiler (with O2 flag for optimization).
Monte Carlo Recipe
I used Monte Carlo to price a Vanilla European Call. It doesn’t have a real sense to perform Monte Carlo on this kind of option but the point of the article is in the execution time.
The procedure is the following:
- Compute the final value of your asset (no need for time discretization since the option isn’t path-dependant, only the final value of the asset is interesting for us here) with a Black Scholes dynamic under risk-neutral probability:
where
- Apply your payoff function to the final value. Here the payoff is .
- Repeat N times.
- Take the mean of all your simulated payoffs and discount the result.
Review of the languages and packages used
Python is a byte-code interpreted language. Your .py code is first compiled into byte-code (.pyc file) and then is interpreted by the Python Virtual Machine.
Cython generates very efficient C code from your Python code only once and then compiles it with C/C++ compilers. It allows to really improve the time-execution of your code without knowing anything about C.
VBA (Virtual Basic for Application): VBA code is compiled into Microsoft P-code (intermediate language) and then executed by a Virtual Machine.
C# is compiled into an intermediate language which is compiled just-in-time (combination of ahead-of-time compilation and interpretation) into the native assembly language of the machine.
C++ is a compiled language. A compiler program translates your code into bytecode that can be processed by the CPU. It is famous for its speed.
Numpy is powerful Python package provides with very efficient N-dimensional arrays and sophisticated functions. It enables to use vectorization instead of loops for some operations which gives a terrible gain in time execution.
Performance assessment
In order to benchmark the performance of the languages, I wrote the same code in each language. I implemented a function Monte_Carlo_Vanilla_Code that runs with 10.000.000 iterations. For each language, I called this function 1, 10, 20, 50, 100 and 200 times and timed the whole execution time.
Code Snippets
Here are the chunks of code that I used in each programming language:
- Python with and without vectorization:
import numpy as np import math import timeit def Black_Scholes(S,R,Vol,T,X): return S*np.exp((R-(Vol**2)/2)*T + Vol*math.sqrt(T)*X) #Uses Numpy Vectorization def Monte_Carlo_Vanilla_Call_Vectorization(S_0,K,R,Vol,T,N_simu): X = np.random.randn(N_simu) #Compute by vectorization: assets_final_value = Black_Scholes(S_0,R,Vol,T,X) payoffs = np.maximum(assets_final_value-K,0) print("Call Price: {}".format(np.mean(payoffs)*math.exp(-R*T))) #Uses Loop def Monte_Carlo_Vanilla_Call(S_0,K,R,Vol,T,N_simu): X = np.random.randn(N_simu) sum_payoffs = 0 for i in range(N_simu): final_value_asset = Black_Scholes(S_0,R,Vol,T,X[i]) sum_payoffs += max(final_value_asset-K,0) print("Call Price: {}".format((sum_payoffs/N_simu)*math.exp(-R*T))) S_0 = 100 K = 100 R = 0.01 T = 1 Vol = 0.2 N_simu = 10000000 start_time = timeit.default_timer() for i in range(1): #Monte_Carlo_Vanilla_Call_Vectorization(S_0,K,R,Vol,T,N_simu) Monte_Carlo_Vanilla_Call(S_0,K,R,Vol,T,N_simu) print(timeit.default_timer() - start_time)
- Cython code (.pyx file). (Observe C variable type declaration). This code is called into another Python file.
import cython cimport numpy as np import numpy as np import math import timeit #Call this one with loop code cdef Black_Scholes(double S,double R,double Vol,double T,double X): return S*np.exp((R-(Vol**2)/2)*T + Vol*math.sqrt(T)*X) #Call this one with vectorization cdef Black_Scholes_Vect(double S,double R,double Vol,double T,np.ndarray[double,mode="c",ndim=1] X): return S*np.exp((R-(Vol**2)/2)*T + Vol*math.sqrt(T)*X) #Uses Numpy Vectorization cpdef Monte_Carlo_Vanilla_Call_Vectorization(double S_0,double K,double R,double Vol,double T,int N_simu): cdef np.ndarray[double,mode="c",ndim=1] X = np.random.randn(N_simu) #Compute by vectorization: cdef np.ndarray[double,mode="c",ndim=1] assets_final_value = Black_Scholes_Vect(S_0,R,Vol,T,X) cdef np.ndarray[double,mode="c",ndim=1] payoffs = np.maximum(assets_final_value-K,0) print("Call Price: {}".format(np.mean(payoffs)*math.exp(-R*T))) #Uses loop cpdef Monte_Carlo_Vanilla_Call(S_0, K, R, Vol, T, N_simu): X = np.random.randn(N_simu) sum_payoffs = 0 final_value_asset = 0 for i in range(N_simu): final_value_asset = Black_Scholes(S_0,R,Vol,T,X[i]) sum_payoffs += max(final_value_asset-K,0) print("Call Price: {}".format((sum_payoffs/N_simu)*math.exp(-R*T)))
- C# code:
static void Main(string[] args) { //Chronometer : var chrono = new Stopwatch(); double vol = 0.3; double R = 0.01; double S0 = 50; double K = 50; double T = 1; //1 year int nb_simulations = 10000000; chrono.Start(); for(int i = 0; i < 10; i++) { Monte_Carlo_Vanilla_Call(S0, T, R, vol, K, nb_simulations); } chrono.Stop(); Console.WriteLine("Time : {0}", chrono.Elapsed); Console.ReadKey(); } static void Monte_Carlo_Vanilla_Call(double S,double T,double R,double vol,double K,int nb_simu) { double final_value_asset; double sum_payoffs = 0; double final_price; //Random number generator: Normal norm = new Normal(0, 1); for (int i = 0; i < nb_simu; i++) { final_value_asset = Generate_Asset(norm.Sample(), S, T, R, vol); sum_payoffs += Math.Max(final_value_asset - K, 0); } final_price = (sum_payoffs/nb_simu) * Math.Exp(-R * T); Console.WriteLine("Price : {0}", final_price); } static double Generate_Asset(double X,double S,double dt,double R, double vol){ return S * Math.Exp((R - (Math.Pow(vol, 2) / 2)) * dt + vol * Math.Sqrt(dt) * X); }
- VBA code:
Option Explicit Sub Monte_Carlo_Vanilla_Call(S_0 As Double, K As Double, R As Double, T As Double, Vol As Double, nb_simulation As Long) Dim sum_payoffs As Long Dim final_asset_value As Double Dim final_price As Double Dim X As Double Dim i As Long sum_payoffs = 0 Randomize For i = 1 To nb_simulation X = Rnd() If X <> 0 Then X = Application.WorksheetFunction.NormInv(X, 0, 1) final_asset_value = Generate_Asset(S_0, R, T, Vol, X) sum_payoffs = sum_payoffs + WorksheetFunction.Max(final_asset_value - K, 0) End If Next i final_price = (sum_payoffs / nb_simulation) * Exp(-R * T) 'Debug.Print (final_price) End Sub Function Generate_Asset(S_0 As Double, R As Double, dt As Double, Vol As Double, X As Double) As Double Generate_Asset = S_0 * Exp((R - (Vol ^ 2) / 2) * dt + Vol * Sqr(dt) * X) End Function Sub main() Dim start As Date Dim S_0 As Double Dim R As Double Dim K As Double Dim T As Double Dim Vol As Double Dim nb_simulation As Long Dim i As Integer R = 0.01 S_0 = 50 K = 50 Vol = 0.3 T = 1 nb_simulation = 1000000 start = Now() For i = 1 To 100 Call Monte_Carlo_Vanilla_Call(S_0, K, R, T, Vol, nb_simulation) Next i Debug.Print (Format(Now() - start, "hh:mm:ss")) 'Debug.Print ("Done") End Sub
- C++ Code (I apologize in advance for experienced C++ programmer, I do not master C++ at all):
#include <iostream> #include <cmath> #include <algorithm> #include <random> #include <chrono> double Generate_Asset(double S,double R,double Vol,double dt,double X){ return S*std::exp((R-(std::pow(Vol,2))/2)*dt + Vol*std::sqrt(dt)*X); } void Monte_Carlo_Vanilla_Call(double S_0,double R,double Vol,double T,double K,double nb_simulation){ static std::default_random_engine e{}; static std::normal_distribution<double> d{0, 1}; double sum_payoffs = 0; double final_asset_value = 0; double X = 0; double call_price; for(int i=0;i<nb_simulation;++i){ X = d(e); final_asset_value = Generate_Asset(S_0,R,Vol,T,X) sum_payoffs += std::max(final_asset_value-K,0.0); } call_price = (sum_payoffs/nb_simulation)*std::exp(-R*T); std::cout << "Price : " << call_price << std::endl; } int main() { double S_0 = 50; double K = 50; double R = 0.01; double Vol = 0.3; double T = 1; double nb_simulation = 10000000; std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now(); for(int i=0;i<100;++i){ Monte_Carlo_Vanilla_Call(S_0,R,Vol,T,K,nb_simulation); } std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now(); std::cout << std::endl << "Execution time : " << std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count() << std::endl; std::cout << "End! Press a key to exit." << std::endl; std::cin.get(); return 0; }
The results
This chart sums up all the results:
Actually, the last line at the bottom of the chart is the overlap of:
- “Cython Numpy Vectorization With Typing”
- “Cython Numpy Vectorization Without Typing”
- “Python Numpy Vectorization”
These last three are therefore the fastest way to perform Monte Carlo here, even faster than C++.
The strength of Numpy:
- Python with loop (orange line): More than 5.000 seconds to perform the 200 calls to the Monte Carlo function.
- Python with Numpy vectorization (brown line): Less than 150 seconds to perform the 200 calls to the Monte Carlo function.
The difference is amazing and shows how powerful is Numpy!
The strength of Cython:
- Python with loop (orange line): More than 5.000 seconds to perform the 200 calls to the Monte Carlo function.
- Cython with loop (green line): This is exactly the same code but generated in C code thanks to Cython. Huge improvement without any change!
- Cython with loop and variable typing (purple line): This is the same code again except that each variable is typed (double, integer, Numpy array…).
Almost 2 times faster just with Cython and variable declaration. - Cython with Numpy vectorization and variables declaration isn’t gaining a lot since Numpy is already really optimized.
C++ (pale green line) is also excellently ranked here and has the expected performance: a little bit less than 300 seconds to perform the 200 calls to the function.
And finally, VBA is slower than C# as expected. Anyway, this is a good performance for VBA, I wasn’t expecting it to be faster than Cython (I wonder if I haven’t done something wrong…).
Conclusion
This benchmark shows the relative performance of several programming languages. I wasn’t aware of the Python’ slowness, it’s pretty disappointing. However, some packages such as Numpy, Cython (and others like Jython) really boost it and make it faster than C++ code!
C++ is undeniably really fast. I understand now why all the books about Monte Carlo pricing deal with C++. However, it’s a language that seems to be pretty hard to master.
VBA and C# give pretty good results but I think it will be interesting to add other programming languages recognized for their swiftness. Therefore, I will try to add soon C, Java, Swift, and Javascript.
By the way, if anyone is interested to provide some results in another language I would be glad to display its code, its result and obviously mention his name.
Thanks, very interesting !
In fact we can see that hard computation tasks such as complex math function are the main things slowing down execution of such algorithms.
By extension, two things are true:
-performance depends on how well native math library is implemented in each languages you’re using in the benchmark,
-perfs are also influenced by the general optimization of the binaries (for compiled languages) by the compilers used, and in the case of java or python it depends on the optimization of the virtual machines that run the byte code (we can see the major optimization between Java’s JVM and Google’s DVM (dalvik virtual machine) which is their own implementation of a java runtime environment for Android operating system).
It might be interesting to see a benchmark of such an algorithm on the Apache Spark environment with Scala.
LikeLiked by 1 person
Thanks a lot Ianis for your comment! Although I am not really aware of Virtual Machine optimization, this is really interesting and think I will try to deepen into this field!
Concerning the math library performance, I believe that it’s more complex than that. The performance will often depend on the data type (and its size) on which you want to apply your mathematical function. An example here: https://stackoverflow.com/questions/3650194/are-numpys-math-functions-faster-than-pythons
LikeLike