Calibration of Heston's stochastic volatility model to an empirical density using a genetic algorithm

Forschung am IVW Köln, 3/2015 Institut für Versicherungswesen Calibration of Heston's stochastic volatility model to an empirical density using a genetic algorithm Urij Dolgov Forschung am IVW Köln, 3 / 2015 Wählen Sie ein Element aus. Dolgov Forschungsstelle FaRis Calibration of Heston's stochastic volatility model to an empirical density using a genetic algorithm Zusammenfassung In diesem Artikel schlagen wir die Verwendung eines genetischen Algorithmus (GA) zur Kalibrierung eines Stochastischen Prozesses an eine empirische Dichte von Aktienrenditen vor. Anhand des Heston Models zeigen wir wie eine solche Kalibrierung durchgeführt werden kann. Neben des Pseudocodes für einen einfachen aber leistungsfähigen GA präsentieren wir zudem auch Kalibrierungsergebnisse für den DAX und den S&P 500. Abstract In this paper we propose the use of genetic algorithms when fitting a stochastic process to the empirical density of stock returns. Using the Heston Model as an example, we show how such a calibration can be carried out. We also present an easy to implement genetic algorithm and provide calibration results for the daily stock returns of the DAX and the S&P 500. Schlagwörter: Aktienrenditen, Dichtefunktion, empirische Dichte, Heston Model, Modellkalibrierung, Stochastische Prozesse Keywords: Empirical Density, Heston Model, Model Calibration, Probability Density, Stochastic Processes, Stock Returns Calibration of Heston’s stochastic volatility model to an empirical density using a genetic algorithm Urij Dolgov 23 December 2014 Abstract In this paper we propose the use of genetic algorithms when fitting a stochastic process to the empirical density of stock returns. Using the Heston Model as an example we show how such a calibration can be carried out. We also present an easy to implement genetic algorithm and provide calibration results for the daily stock returns of the DAX and the S&P 500. 1 The Heston Model and it’s transition density The Heston Model (HM) suggested by Heston (1993) is often seen as the first logical extension of the widely known Black and Scholes (BS) approach. It uses a stochastic volatility instead of the flat one suggested by it’s less sophisticated counterpart. Several empirical studies have shown already that the constant-volatilityassumption contradicts market realities (see e.g. Cont (2001), Guillaume et al. (1997) ) The most salient drawback of the B&S-model is often considered to be it’s inability to replicate the long tails which are observable in daily stockreturns. These however can be captured quite well by Heston’s approach. (see Silva and Yakovenko (2003), Daniel (2003)) The model’s dynamics are characterized by the following three equations dSt=µStdt +√vtStdW(1) t(1) dvt=−γ(vt−θ)dt +κ√vtdW(2) t(2) dW(2) t=ρdW(1) t+p1−ρ2dZt(3) Here Ztis a Wiener process independent of W(1) t.Defining rtby rt= ln(St/S0), applying Ito’s Formula and setting xt=rt−µt one arrives at dxt=−0.5vtdt +√vtdW(1) t(4) 1 A transition density Pt(x, v|vi) for the joint realization of xtand vat time tgiven an initial log-return x= 0 and variance viat t= 0 was constructed by Dragulescu and Yakovenko (2002). However, because the variance is not a directly observable market quantity Pt(x, v|vi) is not suitable for the calibration of the Heston Model to empirical data. Fortunately Dragulescu and Yakovenko (2002) also introduce reduced densities by integrating out the variance. Pt(x|v0) = Z+∞ 0 Pt(x, v|v0)dv (5) This is the density of having log-return xat time tgiven a variance v0at time t= 0. If we wanted to specify a starting variance the function above might be the correct choice. However deciding upon such a value can be quite arbitrary for variance can not be directly observed. To circumvent this uncertainty Dragulescu and Yakovenko (2002) use the stationary density of (2) which is given by Π∗(v) = αα Γ(α) vα−1 θαe−αv/θ, α =2γθ κ2(6) Inserting (6) into (5) and integrating over v0they eventually arrive at the following result: Pt(x) = 1 2πZ+∞ −∞ eipx+Ft(px)dpx(7) with F(t, px) = γθ κ2Γt−2γθ κ2ln cosh Ωt 2+Ω2−Γ2+ 2γΓ 2γΩsinh Ωt 2(8) Γ = γ+ iρκpx(9) Ω = pΓ2+κ2(p2 x−ipx) (10) 2 Empirical density function and fitting Let a series of log-normal stock returns be given by x(∆t) = {x1, . . . , xn}. Here ∆tdenotes the time-step. Thus in the case of 21 trading days per month and ∆t= 1/252, x(1/252) will be a series of daily returns. To make the data compatible with process (4) we also shift every return by µ∆t=∆t nPn i=1 xi and thus reset x(∆t) = {x1−µ∆t, . . . , xn−µ∆t}(11) Furthermore we set xmin = min x(∆t), xmax = max x(∆t) and introduce the bin-size ∆x. The number of bins is then given by M= ceil xmax−xmin ∆x. 2 We then simply determine the relative probabilities pifor each bin by weighting the number of returns in that bin by the total number of observations. For a bin with boundaries [a, b) we also introduce the bin representative ¯xi=a+1 2(b−a). Using the Mrepresentatives {¯x1. . . , ¯xM}we can now define the function pemp(x) = M X i=1 piδ¯xi(x) (12) After the initial empirical function has been constructed one can get rid of outliers by determining the αand 1 −αquantiles. Using these values to reset xmin and xmax one then proceeds to derive a new empirical density (without the outliers). To underline that the shape of the analytical density depends on the choice of (γ, θ, κ, ρ), we set Pt(x) = p(t, x, γ, θ, κ, ρ) The distance function between the empirical distribution and the analytical one is then given by f(γ, θ, κ, ρ) = M X 1=1 [pemp(¯x)−p(∆t, ¯x, γ, θ, κ, ρ))]2(13) To conclude the fitting procedure we have to minimize fwith respect to (γ, θ, κ, ρ). Note that fas defined in (13) is just one choice of a distance function. Alternatives might be the the root-mean-squared error or the absolute distance. 3 The genetic algorithm approach 3.1 Motivation In this section we assume our cost function to be f:Rn→Rand our goal will be to minimize f. The function’s input is thus a vector of the form ~x = (x1, . . . , xn)∈Rnand an input-solution-pair is a tuple (~x, f(~x)) ∈Rn+1. Furthermore we restrict the search-space of our algorithm by limiting it to the hypercube H(~a,~ b) = {~x ∈Rn, ai≤xi≤bi,∀∈{1,...n}} (14) However, if fis set to equal (13) it’s calculation must be carried out via numerical integration and the integral might not be well-behaved for every choice of (γ, θ, κ, ρ). This can lead to the crash of most conventional optimization routines. In the case of the Heston Model fis known explicitly. However, for many stochastic processes the transition density itself is often not known and has to be estimated numerically. Pedersen (1995) and Brandt and Santa-Clara (2001) for example suggest the use of Monte Carlo coupled with a Maximum-Likelihood estimator in oder to construct the transition density. In such a case approaches like Levenberg-Marquard or steepest descent (see Kelley (1999)) would be computationally expensive, for f0will have to be determined numerically. It would be best to use an optimization routine which: 3 •does not depend on the form of f, •doesn’t necessitate the computation of f0, •works reliably for larger n, •is able to disregard local optima •and is inherently parallel which will allow to exploit modern multi-core processors. All of the criteria above are met by genetic algorithms (GA). For an introduction to GA see e.g. Sivanandam and Deepa (2007) or Gen (1997). Despite being considered heuristics they have been repeatedly shown to perform well with complex problems (see e.g. De Jong (1975) , Marco-Blaszka and Desideri (1999)). In the next section we are going to introduce a simple GA containing all the necessary buildings blocks which are inherent to this type of optimization procedure. 3.2 The Algorithm The general structure of most genetic algorithms is quite simple. One first creates an initial population by randomly sampling from the search space. Afterwards the population members are selected and exchange information in order to produce new, fitter solutions. The population is then sorted according to the fitness level and the worst solutions are killed off. The members of this new population are again selected for ”mating” and the entire procedure is repeated. To embed our optimization problem into the GA-framework we must introduce some definitions first. From here on ~x will be called a chromosome and the tuple (~x, f(~x)) ∈Rn+1 a member of the population ˜ Pwhich is simply a set of such tuples. To make the notation more compact we will refer to a population member via ˜m, it’s chromosome via ˜m(1) and it’s cost via ˜m(2) . In this context ˜ Piwill be the i-th population member. Here we also access the chromosome via ˜ Pi(1) ∈Rnand the cost via ˜ Pi(2) ∈R. To initialize a GA one first has to generate an initial population by randomly sampling Nini points from the space H(~a,~ b) (see Algorithm 1) 4 input : The dimension of the objective function’s input n, the size of the initial population Nini, an empty container ˜ Pini and the vectors ~a,~ b∈Rn output: The container ˜ Pini filled with Nini elements local :~x ∈Rn for i= 1 to Nini do for i= 1 to ndo xi=drawUniformRandom [ai,bi]; end Add (~x, f(~x)) to ˜ Pini; end Algorithm 1: Creating an initial population The next step consists of deciding which members of the population should mate. Here the best solution should also have the best chances of passing on it’s information. To achieve that we use a cost-weighted-selection approach. Thus the lower the cost of a given population member the higher it’s probability to be selected for mating. Assuming that the current population ˜ Pis already sorted so that ˜ P1(2) ≤˜ Pi(2) ∀˜ Pi∈˜ Pthe selection probabilities are given by pi= Pi(2) −PN(2) PN j=1(Pj(2) −PN(2))∀i∈ {1, . . . , n}(15) Using this probabilities one than constructs a discrete probability distribution with P(X=i) = piand uses it to draw members from ˜ P. The number of mating-pair selections is regulated by the variable η∈[0,1] specifying the fraction of the population that will mate (see Algorithm 2) Note that algorithm 2 allows self-pairing. Thus a member can be paired with itself. 5 input : The current population ˜ Pand it’s size Nand η∈R output: A list of mating-pairs ˜ M local :~p ∈RN, ˜m1∈Rn+1 , ˜m2∈Rn+1 ˜ P=sortPopulation (˜ P) ; // so that ˜ P1(2) ≤˜ Pi(2) ∀˜ Pi∈˜ P for i= 1 to Ndo pi= Pi(2)−PN(2) PN j=1 (Pj(2)−PN(2))  end for i= 1 to b0.5Nηcdo ˜m1=drawMemberFromPopulation (~p); ˜m2=drawMemberFromPopulation (~p); Add ( ˜m1,˜m2) to ˜ M; end Algorithm 2: Selecting members for mating Now that we have selected the mating pairs, we proceed to introduce another crucial genetic operator - the crossover. The crossover constitutes the step where the exchange of information between the different solutions present in the population takes place. The key idea is to combine data of two members (to ”mate” them) in order to produce two fitter members. How this is accomplished varies depending on the concrete problem. The approach used here is described in Algorithm 3. Using the mating pairs obtained via Algorithm 2 one goes on to combine the genetic data of the two members to produce two new ones. The new members inherit most of their parent’s data but one randomly selected gene. We collect this results in a container ˜ Pch and set ˜ P=˜ P∪˜ Pch. To avoid an exponential population growth we introduce the variable Nmax, sort ˜ Pso that ˜ P1(2) ≤˜ Pi(2) ∀˜ Pi∈˜ Pand remove all ˜ Piwith i>Nmax. This way the population size is kept constant to reduce computational overhead. Like many optimization routines the GA-approach might also suffer from premature convergence. Thus it finds one local minimum and converges to it neglecting the better optimal solutions within the search-space. In the context of the GA this is often caused by the genetic-drift. This happens when the population ˜ Pcontains one member ˜mbest which is significantly fitter than the others. ˜mbest will thus be frequently selected for mating and will end up dominating the ”gene-pool”. After several iterations the entire population might end up converging to ˜mbest . The usual approach to remedy this situation is to introduce a steady flow of new and unbiased genetic information to the population. This can be achieved by randomly sampling points from the search space and adding them to ˜ P. Thus one runs Algorithm 1 to produce a specified number of new members Nmut and adds them to the current population. This 6 is also often referred to as introducing mutation. Here it is important to note that this new members will be added to ˜ Pafter it has been truncated to Nmax. This way they will have a chance to mate in the next run. input : A set of mating pairs ˜ M, their number NM, factor β∈[0,1], the current population ˜ Pand Nmax output: An updated population ˜ P local : ˜m1∈Rn+1 , ˜m2∈Rn+1 ,a∈N,~x, ~y ∈Rn,z1, z2∈R, ~z1, ~z2∈Rn,Nnew for i= 1 to NMdo ( ˜m1,˜m2) = ˜ Mi; a=drawInteger ({1, . . . , n}); // Uniform random draw of an integer from {1, . . . , n} ~x = ˜m1(1) = {x1, . . . , xa,...xn}; ~y = ˜m2(1) = {y1, . . . , ya,...yn}; z1=xa−β(xa−ya); z2=ya+β(xa−ya); ~z1={x1, . . . , xa−1, z1, xa+1 . . . xn}; ~z2={y1, . . . , ya−1, z1, ya+1 . . . yn}; ˜m1= (~z1, f(~z1)) ˜m2= (~z2, f(~z2)) Add ˜m1and ˜m2to ˜ Pch; end ˜ P=˜ P∪˜ Pch; ˜ P=sortPopulation (˜ P) ; // so that ˜ P1(2) ≤˜ Pi(2) ∀˜ Pi∈˜ P Nnew=getPopSize (˜ P); ˜ P=˜ P\{˜ PNmax+1,..., ˜ PNnew }; Algorithm 3: Mating/crossover and regulating population size Finally we have to define a stopping-criteria for the algorithm. The most straight-forward approach and the one we use here, is to stop when the fittest population member ˜ P1hasn’t changed K-times in a row. 7 References Michael W. Brandt and Pedro Santa-Clara. Simulated likelihood estimation of diffusions with an application to exchange rate dynamics in incomplete markets. NBER Technical Working Paper 0274, National Bureau of Economic Research, Inc, 2001. URL https://ideas.repec.org/p/nbr/nberte/0274.html. Rama Cont. Empirical properties of asset returns: stylized facts and statistical issues. Quantitative Finance, 1:223–236, 2001. Gilles Daniel. Goodness-of-fit of the heston model. arXiv:cs/0305055, May 2003. URL http://arxiv.org/abs/cs/0305055. Kenneth Alan De Jong. An Analysis of the Behavior of a Class of Genetic Adaptive Systems. PhD thesis, University of Michigan, Ann Arbor, MI, USA, 1975. AAI7609381. Adrian A. Dragulescu and Victor M. Yakovenko. Probability distribution of returns in the heston model with stochastic volatility. arXiv:cond-mat/0203046, March 2002. URL http://arxiv.org/abs/cond-mat/0203046. Quantitative Finance 2, 443 (2002). Mitsuo Gen. Genetic Algorithms and Engineering Design. John Wiley & Sons, January 1997. ISBN 9780471127413. Dominique M. Guillaume, Michel M. Dacorogna, Rakhal R. Dav, Ulrich A. Mller, Richard B. Olsen, and Olivier V. Pictet. From the bird’s eye to the microscope: A survey of new stylized facts of the intra-daily foreign exchange markets. Finance and Stochastics, 1(2):95– 129, April 1997. ISSN 0949-2984. doi: 10.1007/s007800050018. URL http://link.springer.com/article/10.1007/s007800050018. S. L. Heston. A closed-form solution for options with stochastic volatility with applications to bond and currency options. Review of Financial Studies, 6(2): 327–343, January 1993. ISSN 0893-9454, 1465-7368. doi: 10.1093/rfs/6.2.327. URL http://rfs.oxfordjournals.org/content/6/2/327. C. T. Kelley. Iterative Methods for Optimization. SIAM, January 1999. ISBN 9780898714333. Nathalie Marco-Blaszka and Jean-Antoine Desideri. Numerical solution of optimization test-cases by genetic algorithms. February 1999. URL http://hal.inria.fr/inria-00073055. Morten Bjerregaard Pedersen. A new approach to maximum likelihood estimation for stochastic differential equations based on discrete observations. Scandinavian Journal of Statistics, 22:55–71, 1995. 14 A. Christian Silva and Victor M. Yakovenko. Comparison between the probability distribution of returns in the heston model and empirical data for stock indexes. Physica A: Statistical Mechanics and its Applications, 324(1): 303–310, 2003. ISSN 0378-4371. S. N. Sivanandam and S. N. Deepa. Introduction to Genetic Algorithms. Springer Science & Business Media, October 2007. ISBN 9783540731900. 15 Impressum Diese Veröffentlichung erscheint im Rahmen der Online-Publikationsreihe „Forschung am IVW Köln“. Alle Veröffentlichungen dieser Reihe können unter www.ivw-koeln.de oder hier abgerufen werden. Forschung am IVW Köln, 3/2015 Dolgov: Calibration of Heston's stochastic volatility model to an empirical density using a genetic algorithm Köln, Februar 2015 ISSN (online) 2192-8479 Herausgeber der Schriftenreihe / Series Editorship: Prof. Dr. Lutz Reimers-Rawcliffe Prof. Dr. Peter Schimikowski Prof. Dr. Jürgen Strobel Institut für Versicherungswesen / Institute for Insurance Studies Fakultät für Wirtschaftsund Rechtswissenschaften / Faculty of Business, Economics and Law Fachhochschule Köln / Cologne University of Applied Sciences Web www.ivw-koeln.de Schriftleitung / Contact editor’s office: Prof. Dr. Jürgen Strobel Tel. +49 221 8275-3270 Fax +49 221 8275-3277 Mail [email protected] Institut für Versicherungswesen / Institute for Insurance Studies Fakultät für Wirtschaftsund Rechtswissenschaften / Faculty of Business, Economics and Law Fachhochschule Köln / Cologne University of Applied Sciences Gustav Heinemann-Ufer 54 50968 Köln Kontakt Autor / Contact author: Urij Dolgov Institut für Versicherungswesen / Institute for Insurance Studies Fakultät für Wirtschaftsund Rechtswissenschaften / Faculty of Business, Economics and Law Fachhochschule Köln / Cologne University of Applied Sciences Gustav Heinemann-Ufer 54 50968 Köln Tel. +49 221 8275-3271 Fax +49 221 8275-3277 Mail [email protected]e Zuletzt erschienen im Rahmen von „Forschung am IVW Köln“ 2015  Heep-Altiner, Berg: Mikroökonomisches Produktionsmodell für Versicherungen, Nr. 2/2015  Institut für Versicherungswesen: Forschungsbericht für das Jahr 2014, Nr. 1/2015 2014  Müller-Peters, Völler (beide Hrsg.): Innovation in der Versicherungswirtschaft, Nr. 10/2014  Knobloch: Zahlungsströme mit zinsunabhängigem Barwert, Nr. 9/2014  Heep-Altiner, Münchow, Scuzzarello: Ausgleichsrechnungen mit Gauß Markow Modellen am Beispiel eines fiktiven Stornobestandes, Nr. 8/2014  Grundhöfer, Röttger, Scherer: Wozu noch Papier? Einstellungen von Studierenden zu EBooks, Nr. 7/2014  Heep-Altiner, Berg (beide Hrsg.): Katastrophenmodellierung - Naturkatastrophen, Man Made Risiken, Epidemien und mehr. Proceedings zum 6. FaRis & DAV Symposium am 13.06.2014 in Köln, Nr. 6/2014  Goecke (Hrsg.): Modell und Wirklichkeit. Proceedings zum 5. FaRis & DAV Symposium am 6. Dezember 2013 in Köln, Nr. 5/2014  Heep-Altiner, Hoos, Krahforst: Fair Value Bewertung von zedierten Reserven, Nr. 4/2014  Heep-Altiner, Hoos: Vereinfachter Nat Cat Modellierungsansatz zur Rückversicherungsoptimierung, Nr. 3/2014  Zimmermann: Frauen im Versicherungsvertrieb. Was sagen die Privatkunden dazu?, Nr. 2/2014  Institut für Versicherungswesen: Forschungsbericht für das Jahr 2013, Nr. 1/2014 2013  Heep-Altiner: Verlustabsorbierung durch latente Steuern nach Solvency II in der Schadenversicherung, Nr. 11/2013  Müller-Peters: Kundenverhalten im Umbruch? Neue Informationsund Abschlusswege in der Kfz-Versicherung, Nr. 10/2013  Knobloch: Risikomanagement in der betrieblichen Altersversorgung. Proceedings zum 4. FaRis & DAV-Symposium am 14. Juni 2013, Nr. 9/2013  Strobel (Hrsg.): Rechnungsgrundlagen und Prämien in der Personenund Schadenversicherung - Aktuelle Ansätze, Möglichkeiten und Grenzen. Proceedings zum 3. FaRis & DAV Symposium am 7. Dezember 2012, Nr. 8/2013  Goecke: Sparprozesse mit kollektivem Risikoausgleich - Backtesting, Nr. 7/2013  Knobloch: Konstruktion einer unterjährlichen Markov-Kette aus einer jährlichen MarkovKette, Nr. 6/2013  Heep-Altiner et al. (Hrsg.): Value-Based-Management in Non-Life Insurance, Nr. 5/2013  Heep-Altiner: Vereinfachtes Formelwerk für den MCEV ohne Renewals in der Schadenversicherung, Nr. 4/2013  Müller-Peters: Der vernetzte Autofahrer – Akzeptanz und Akzeptanzgrenzen von eCall, Werkstattvernetzung und Mehrwertdiensten im Automobilbereich, Nr. 3/2013  Maier, Schimikowski: Proceedings zum 6. Diskussionsforum Versicherungsrecht am 25. September 2012 an der FH Köln, Nr. 2/2013  Institut für Versicherungswesen: Forschungsbericht für das Jahr 2012, Nr. 1/2013