Computational techniques in queueing and fluctuation theory
Mohammad Asghari, N.
Publication date
2014
Document Version
Final published version
Link to publication
Citation for published version (APA):
Mohammad Asghari, N. (2014). Computational techniques in queueing and fluctuation theory.
General rights
It is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s) and/or copyright holder(s), other than for strictly personal, individual use, unless the work is under an open content license (like Creative Commons).
Disclaimer/Complaints regulations
If you believe that digital publication of certain material infringes any of your rights or (privacy) interests, please let the Library know, stating your reasons. In case of a legitimate complaint, the Library will make the material inaccessible and/or remove it from the website. Please Ask the Library: https://uba.uva.nl/en/contact, or a letter to: Library of the University of Amsterdam, Secretariat, Singel 425, 1012 WP Amsterdam, The Netherlands. You will be contacted as soon as possible.
Computational techniques
in queueing and fluctuation theory
Naser Mohammad Asghari
Naser
Mohammad
Asghari
for the public defense
of my PhD thesis
Computational techniques
in queueing and fluctuation
theory
that will take place on
Tuesday
November 25, 2014
at 14.00
in the Agnietenkapel,
Oudezijds voorburgwal 231,
Amsterdam
The defense will be
followed by a reception
Computational techniques
in queueing and fluctuation theory
Computational techniques
in queueing and fluctuation theory Naser Mohammad Asghari
Korteweg-de Vries Instituut voor Wiskunde
Faculteit der Natuurwetenschappen, Wiskunde en Informatica Proefschrift Universiteit van Amsterdam
Copyright c 2014 by Naser Mohammad Asghari, Amsterdam All rights reserved. No part of this book may be reproduced, in any form or by any means, without permission in writing from the author.
Computational techniques
in queueing and fluctuation theory
ACADEMISCH PROEFSCHRIFT
ter verkrijging van de graad van doctor aan de Universiteit van Amsterdam
op gezag van de Rector Magnificus prof. dr. D.C. van den Boom
ten overstaan van een door het college voor promoties ingestelde commissie,
in het openbaar te verdedigen in de Agnietenkapel op dinsdag 25 november 2014, te 14:00 uur
door
Naser Mohammad Asghari
Promotiecommissie
Promotor: Prof. dr. M. R. H. Mandjes Overige leden: Dr. P. J. C. Spreij
Prof. dr. R. Núñez Queija Prof. dr. ir. C. W. Oosterlee Prof. dr. J. H. van Zanten Prof. dr. D. T. Crommelin Prof. dr. B. F. Heidergott
Acknowledgments
When I was school student, I was fascinated by almost every topic, mathematics, physics, chemistry, biology, history, geography, ... (except arts!), and I wanted to be expert in all of them. One day I studied physics, the other day history, and so on. Finally I found out I like physics and mathematics most. So I started mathematics at university, but the next year I changed to physics. After I got graduated in physics, I was introduced to economics and finance by my brother Mohsen, and my best friend, Mehdi. Physicists had already been active in economy and finance for some time, and they called their research econophysics. Researching in econophysics led me to stochastic processes and mathematical finance. I was very lucky that I managed to persuade Michel Mandjes to supervise me, despite the fact that I did not have a mathematics background. Although I had studied probability and stochastic processes, my knowledge was not enough to start doing research. Michel trusted me and guided me such that I managed to get on the right track in quite a short period of time. In fact, it would not be possible to accomplish this thesis without his great help and support. Michel, thank you for your trust, support, guidance and kindness. You also supported me apart from my thesis and I always appreciate it.
I would also like to thank Peter Spreij. He gave me the chance to be an assistant for his courses at UvA. He is a great teacher and those courses still help me in my research and in my job. He also supported me when I was looking for a job. Thank you, Peter.
I met Martijn Pistorius during a summer school in mathematical finance, in 2011 in Ljubljana. I recall that I enjoyed his lectures. When he came to KdVI, whenever I had challenges in my research, he always had valuable suggestions and comments. I would like to thank him. KdVI is a prestigious institute, and I am really proud of being a member of it. I would not have gotten this chance without the support of Jan Wiegerinck, to whom I am very grateful.
I would also like to thank Evelien, Henneke and Marieke for their help. I never felt lonely at KdVI with Nabi, Paul, Jevgenijs, Ricardo, Enno, Loek, Arie, and Piotr. I will never forget the dinners we had at Jevgenijs’s place.
I am also so thankful to Peter den Iseger, Anwar Walid and Krzysztof Debicki for the nice collaboration when writing our joint papers.
I would like to thank Drona Kandhai for giving me the chance to work in his team at ING. I really enjoy working with my managers and colleagues at ING: Dirk, Veronica, Bart, Geert-Jan, Joanna, Artem, Dmytro, Daan, Jef, Markus, Xiaoyu, Carlos and Frits Koen.
Outside the University and ING, I belong to an Iranian community. We have a great time together. I would like to thank Mohammad & Maryam, Vahid & Sara, Hodjat & Marzieh, Mehdi & Marzieh, Afshin & Fatemeh, Ali & Azadeh, Mahdi & Mahdieh, Masoud & Hoda, Amin & Fatemeh, Naser & Aylar, Shayan & Narges, Mohammad & Fahimeh, Saeed & An-disheh, Farzin & Mahsa, Hossein, Mahdi, Jafar, Narges, Rojman, Mohammad, Abbas, Da-nial, Amir, Behnaz, Maryam and Rokhsareh.
I want to thank my family. How would I get here without your love, support and encour-agement? My mother and father, sisters and brother, nephew and nieces, my father, mother, brothers and sisters in law. I love you all and I cannot imagine a moment without you. My wife, Sareh, deserves a special acknowledgment. Thank you for your love and support. You always give me energy and encouragement to continue. This thesis is dedicated to you.
Naser Mohammad Asghari, Amsterdam, October 2014
Contents
Acknowledgments v List of Figures ix List of Tables xi 1 Introduction 1 1.1 Motivation . . . 1 1.2 Models . . . 31.3 Preliminaries on Lévy fluctuation theory . . . 4
1.4 Preliminaries on Markov fluid models . . . 15
1.5 Outline of thesis, contributions . . . 20
2 Numerical techniques in Lévy fluctuation theory 25 2.1 Introduction . . . 25
2.2 Preliminaries . . . 29
2.3 Laplace inversion . . . 30
2.4 Approximation with rational Laplace transform . . . 36
2.5 Small jumps . . . 38
2.6 Beta processes . . . 40
2.7 Discussion and concluding remarks . . . 43
3 Evaluation of option prices 51 3.1 Introduction . . . 51
3.2 Preliminaries . . . 54
3.3 Transforms of prices and Greeks of lookback options . . . 61
3.4 Numerical validation . . . 65
3.5 Discussion . . . 78
3.6 Concluding Remarks . . . 80
4 Asymptotics of the supremum of a Lévy process 81 4.1 Introduction . . . 81
4.2 Asymptotics . . . 82
4.3 Importance sampling . . . 84
5 Energy-Efficient Scheduling in Multi-Core Servers 89 5.1 Introduction . . . 89
5.2 Energy Cost Function and Models . . . 92
5.3 Optimization of Energy Consumption . . . 100
5.4 Robustness Analysis . . . 108
5.5 Conclusion . . . 112
Bibliography 115
Samenvatting 122
Summary 124
About the author 126
List of Figures
5.1 Optimal cost of static serving and sleep mode strategies vs. buffering cost rate. 102 5.2 Thresholds in the 1-threshold and the hysteretic strategy with respect to β . . 103 5.3 Optimal cost of the hysteretic strategy and the 1-threshold strategy. . . 103 5.4 Optimal cost of strategies with many thresholds. . . 106 5.5 Service rate as function of buffer content in continuous strategy. . . 106
List of Tables
2.1 Brownian Motion. . . 34
2.2 Compound Poisson with exponential jumps. . . 46
2.3 Compound Poisson with Weibull jumps. . . 46
2.4 Compound Poisson with Weibull jumps. . . 47
2.5 Compound Poisson with Pareto jumps. . . 47
2.6 Compound Poisson with shifted-Pareto jumps. . . 48
2.7 Compound Poisson with both upward and downward jumps. . . 48
2.8 CGMY-like upward-jumps. . . 49
2.9 Variance Gamma process. . . 49
2.10 Beta process. . . 50
3.1 Black-Scholes model. . . 68
3.2 Greeks corresponding to Black-Scholes model. . . 69
3.3 Merton model. . . 70
3.4 Kou model. . . 72
3.5 CGMY model. . . 75
3.6 Greeks corresponding to CGMY model. . . 76
3.7 Beta model. . . 77
4.1 Simulation results corresponding to Compound Poisson process. . . 86
4.2 Simulation results corresponding to Variance Gamma process. . . 87
5.1 Optimal cost of single-server strategies. . . 104 xi
5.2 Optimal cost of multiple-server strategies. . . 107 5.3 Robustness with respect to changes in the mean on-time. . . 109 5.4 Robustness with respect to changes in the distribution of the on-times, with
the mean on-time unchanged. . . 110 5.5 Robustness of multiple-server strategies with respect to changes in the
distri-bution of the on-times, with the mean on-time unchanged. . . 111 5.6 Robustness with respect to changes in the distribution of the on-times, with
the mean on-time changing as well. . . 111 5.7 Robustness of multiple-server strategies with respect to changes in the
distri-bution of the on-times, with the mean on-time changing as well. . . 112
Chapter
1
Introduction
1.1
Motivation
Stochastic phenomena are everywhere around us. Consider the example of the valuation of options (or other financial products) in financial markets. The future prices of the un-derlying assets being highly uncertain, we need a model to describe these. Such a model should be able to accurately capture the random features of the underlying stochastic pro-cess. Realizing that option prices are essentially functionals of the evolution of the price of the underlying asset, the model can be used to price the options.
Another example in which randomness plays an important role, is in the design of vari-ous types of communication infrastructures. A commonly used paradigm is to model such systems as queueing networks. The design objective typically reflects the tradeoff between the cost and the quality-of-service delivered to the customers; adding capacity improves the perceived quality, but obviously comes at a price. In recent years substantial emphasis has been put on designing the network such that it efficiently uses energy resources, for instance by adaptively changing the processor speeds of the queues, as a function of the current workload. Such algorithms need to be set up delicately, as they should not compromise the performance perceived by the network’s users.
These are just illustrative examples of instances where stochastic modeling offers a useful mathematical framework. Of course in many other situations such models can be applied as well. One could think of a wide variety of other domains, such as population biology, chemical reaction networks, and statistical physics. The use of techniques from stochastic
2 1.1. MOTIVATION modeling, to optimize the efficient usage of network resources, obviously extends beyond the application of communication networks that we mentioned above: also in transport net-works and logistic netnet-works this approach is intensively used. The list of potential appli-cation areas keeps on expanding; recently developed areas include social networks, brain modeling, and forensic sciences.
Modeling real-life situations by means of stochastic systems is typically a first step; in or-der for these models to have practical use, they need to allow fast and accurate (numerical) evaluation. To illustrate this, let us go back to the two motivating examples that we intro-duced above. In the context of option pricing, traders need to respond to the clients’ requests nearly instantly, and this requires that they need to be able to compute prices virtually in real time. The alternative is to rely on simulation-based techniques, which are fast only when specific simplifications have been imposed on the underlying model (for instance by assum-ing a simplistic volatility model, or by assumassum-ing the underlyassum-ing stochastic process is just Brownian); when the ambition is to use more sophisticated models, simulation techniques become inherently problematic. This motivates the research on fast and accurate numerical computation techniques that do not require us to a priori simplify the underlying stochastic processes.
Also in the setting of the (optimal) design of energy-efficient communication networks, com-putational issues play a prominent role. A typical objective function includes the energy consumption per time unit as well as a performance measure that reflects the quality-of-service, and the idea is to find the service policy that strikes an optimal balance between these components (or the less ambitious objective of identifying a policy that is ‘close’ to this optimum). The optimization tries to find, for any specific condition the network can be in, what the best service strategy is. The complexity of finding an optimum lies in the huge parameter space that we have to optimize over, as well as the sometimes unpleasant specific properties of the objective function (it is not necessarily a smooth unimodal function, for instance). In the setting of designing a network the optimization routine can in principle be performed off-line. It should be realized, however, that the underlying algorithm requires that the objective function be evaluated potentially very frequently, and as a consequence we need a technique that evaluates the system’s performance (for a single parameter instance) fast and sufficiently accurately. Perhaps equally important is that it is made sure that the resulting design is robust: when the input parameters differ slightly from their estimated values, the system should still behave ‘nearly’ optimally.
CHAPTER 1. INTRODUCTION 3
1.2
Models
A canonical model in stochastic modeling is the so-called random walk. Consider a sequence of independent and identically distributed increments Y1, Y2, . . ., and the associated random
walk Sn:= n i=1 Yi.
The probabilistic behavior of these random walks has been an object of intensive study. Classical results are the law of large numbers, saying that Sn/n converges (in various
spe-cific forms) to the corresponding mean, and the central limit theorem, stating that a centered and normalized version of Sn/nconverges to a standard Normal random variable.
In the context of the option pricing example introduced above, observations from financial data suggest that models with independent and identically distributed increments form a natural framework. On the other hand, as in this example time is in principle not slotted, it is less appropriate to consider a discrete-time framework. This motivates why it makes sense to look at the continuous-time counterpart of the above random-walk model, viz. the Lévy model. A Lévy process (Xt)t≥0 is a real-valued continuous-time process, with X0 = 0,
such that all increments are independent and identically distributed (meaning that for all s, hand t we have that Xs+h−Xsand Xt+h−Xthave the same distribution, and that Xt+h−Xt
is independent of Xt). In some cases processes in finance are not well modeled by (Xt)t≥0
itself, but rather by (eXt)
t≥0 (thus also avoiding the sometimes problematic issue that the
underlying process can attain negative values).
In the example of option pricing, various payoff structures involve both the value XT of the
underlying at maturity T and the largest value ¯XT attained until T . An example here is
the so-called barrier option: there, for a given strike K and barrier H, the payoff is given by max{XT − K, 0} provided that ¯XT ≥ H; if we have insight into the probabilistic properties
of these objects, we are in a position to price the option. This explains the interest in getting a handle on the joint distribution of XTand the running maximum ¯XT. As it turns out, for Lévy
processes there is a vast body of literature available on this joint distribution, commonly referred to as Wiener-Hopf theory. We recapitulate the first principles, as well as a collection of key results, of this theory in Section 1.3.
4 1.3. PRELIMINARIES ON LÉVY FLUCTUATION THEORY cannot be processed immediately upon arrival is temporarily stored in a buffer, thus nat-urally leading to a notion of a queueing process. The Lévy process defined above can in principle attain any real value (i.e., positive and negative), so to make it a queueing process we should prevent it from becoming negative. The (commonly followed) way of truncating the process is by applying a so-called reflection or regulation: we introduce a queueing process (or: workload process) by
Qt:= sup s≤t
(Xt− Xs);
observe that this entails that the workload is a functional of the driving Lévy process (Xt)t≥0.
It is remarked that in such a Lévy-driven queue the amount of traffic fed into the system in two disjoint periods of time is independent, which is typically not true in the context of communication networks: if a user is generating traffic at a high rate at some point in time, it is relatively likely he or she is still doing that some small amount of time later. This observation motivates why traffic is often modeled as a so-called on-off process: the users activity level alternates between generating traffic at some rate r for a while, and being silent, thus creating some positive dependence within the queue’s input process. Subsequent on-and off-times may be assumed to be independent.
1.3
Preliminaries on Lévy fluctuation theory
Brownian motion is one of the most frequently used stochastic models, being applied in a va-riety of domains. Among its most important properties are the continuity of its sample paths and its scale invariance. It is noted, however, that many phenomena which are modeled by Brownian motion, do not have these properties. In finance, for example, returns are often modeled by Brownian motion, but they tend to jump up and down; these discontinuities can not be ignored in many situations [31].
Another classical stochastic process is the Poisson process. This has discontinuous sample paths (as it jumps up by steps of size 1, after exponentially distributed times). A Poisson process is a non-decreasing process, and therefore has paths of bounded variation over finite time (as opposed to Brownian motion).
Although Brownian motion and the Poisson process are seemingly completely different, they share a number of important features. Both have right-continuous paths with left limits, and they both are Markovian processes that start at the origin. They belong to a wide class of stochastic processes, Lévy processes, named after the French mathematician Paul Lévy, which
CHAPTER 1. INTRODUCTION 5 play an important role in this thesis. We primarily focus on these processes with a specific focus on the analysis of the probabilistic properties related to its extreme values; this branch of research is commonly referred to as fluctuation theory. For a complete treatment of Lévy processes and fluctuation theory, we refer to e.g. the textbook [68].
A process X ={Xt, t ≥ 0} defined on a probability space (Ω, F, P) is a Lévy process if it
has right-continuous paths with left limits, and if it has stationary and independent incre-ments, with X0 = 0. ‘Stationarity’ in this context means that increments corresponding to
a fixed time interval are identically distributed; ‘independence’ refers to the property that increments corresponding to non-overlapping time intervals behave statistically indepen-dently.
An important, and highly convenient, property of the Lévy process is that its characteristic function obeys a closed-form formula, which is in fact an immediate consequence of the pro-cess having stationary and independent increments. From the definition of a Lévy propro-cess, it is known that Xtis infinitely divisible for any t > 0. Realizing that, for any n ∈ N, Xt
equals
Xt= Xt/n+ (X2t/n− Xt/n) +· · · (Xt− X(n−1)t/n),
it follows that
logEeisXt = t logEeisX1
= tξ(s)
(first for rational s, and with a limiting argument also for real s). The function ξ(s) := logEeisX1 is often referred to as the characteristic exponent of the Lévy process. As stated
in the following theorem, it can be shown that a Lévy process can be uniquely defined by a triple (μ, σ, Π) [68, 34, 21].
Theorem 1.3.1. Lévy-Khintchine formula for Lévy processes. Suppose that μ ∈ R, σ ≥ 0. Let the measure Π be concentrated onR\{0}, in such a way that the regularity condition
Rmin{1, x2}Π(dx) < ∞ is met. This triple (μ, σ, Π) defines for any s ∈ R,
ξ(s) = logEeisX1 = iμs−1
2s 2 σ2+ ∞ −∞(e isx− 1 − isx1 {|x|<1})Π(dx).
Then there exists a probability space (Ω, F, P) on which a Lévy process can be defined hav-ing Lévy characteristic exponent ξ(s). Π is called the correspondhav-ing Lévy measure, while μ is often referred to as the drift, and σ2
6 1.3. PRELIMINARIES ON LÉVY FLUCTUATION THEORY The characteristic exponent can be rewritten as
ξ(s) = iμs−1 2s 2 σ2 + Π(R\(−1, 1)) |x|≥1(e isx− 1) Π(dx) Π(R\(−1, 1)) + 0<|x|<1 (eisx− 1 − isx)Π(dx) .
It should be mentioned that in case Π(R\(−1, 1)) = 0 the second term is left out. Based on this formula of the characteristic exponent, the Lévy process Xtcan be decomposed as the
independent sum of processes X(1)
, X(2)
and X(3)
which are described as follows (Lévy-Itô decomposition). In the first place, X(1)
is a linear Brownian motion with drift. Then, X(2)
is a compound Poisson process with rate Π(R\(−1, 1)), where the jumps are independent and identically distributed with distribution Π(dx)/Π(R\(−1, 1)) concentrated on {x : |x| ≥ 1}. Finally, concentrating on the last term, it is first observed that it can be written as
0<|x|<1 (eisx− 1 − isx)Π(dx) = n≤0 λn 2−(n+1)≤|x|<2−n (eisx− 1)Fn(dx)− isλn 2−(n+1)≤|x|<2−n xFn(dx)
where λn := Π(2−(n+1) ≤ |x| < 2−n)and Fn(dx) := Π(dx)/λn. The component X(3) can
thus be considered as the superposition of (at most) a countable number of independent compound Poisson processes with different rates and linear drift. In fact, X(3)
is a square integrable martingale with an almost surely countable number of jumps on each finite time interval. Importantly, the number of jumps can be infinite almost surely, leading to the class of Lévy models with infinite activity; we will intensively work with this class later in this thesis.
Further relevant classifications are the following. If Π(−∞, 0) = 0, then it follows from the Lévy-Itô decomposition that the corresponding Lévy process has no negative jumps. In this case it is referred as a spectrally positive Lévy process. On the contrary, a Lévy process is called spectrally negative if−X is spectrally positive (i.e., it has no positive jumps). These two classes of processes are generally indicated byS+ andS− respectively, and referred to as
CHAPTER 1. INTRODUCTION 7 one-sided Lévy processes often very explicit analysis is possible.
Let X be a spectrally positive process, and assume in addition(0,∞)max(1, x)Π(dx) <∞,
σ = 0and μ≥ 0. Then, again from the Lévy-Itô decomposition, it follows that the process X has non-decreasing paths. Such a process is referred to as a subordinator. Given a Lévy process Xtand an independent subordinator τs, we can introduce another Lévy process by
sampling Xtat stochastic time epochs which are defined by the subordinator. More precisely,
suppose that Xtis a Lévy process with characteristic exponent ξ and τ = {τs : s ≥ 0} is
an independent subordinator with characteristic exponent Ξ. Then the process Y , which is defined by Xτs, is a Lévy process, and its characteristic exponent is given by Ξ◦ ξ [68]. For
example, a possible representation of the so-called Variance Gamma process (a frequently used infinite-activity Lévy process) corresponds to sampling a Brownian motion at times that result from a Gamma process [34].
1.3.1
Wiener-Hopf factorization
A collection of important results in Lévy fluctuation theory are immediate consequences of the so-called Wiener-Hopf factorization. In fact, this Wiener-Hopf factorization provides a powerful tool with several applications in probability (e.g. related to finance). In this sec-tion we first introduce a few relevant concepts, and then we roughly sketch a proof for the Wiener-Hopf factorization theorem in a discrete-time framework. The continuous-time set-ting is considerably more technical, and therefore we decided to just state the main result, and leave the underlying considerations out.
As it was mentioned, any Lévy process is defined on a probability space (Ω, F, P). Let F be the filtrationF = {Ft : t≥ 0}, so that we obtain a filtered probability space (Ω, F, F, P), on
which we assume X is defined. Then the non-negative random variable τ , defined on the same filtered probability space, is called stopping time if{τ ≤ t} ∈ Ftfor all t > 0. It should
be mentioned that it is not a priori ruled out that a stopping time could have the property thatP(τ = ∞) > 0. Now suppose that τ is a stopping time. The process ˜X ={ ˜Xt : t≥ 0}
where
˜
Xt= Xτ +t− Xτ, t≥ 0
defined on{τ < ∞} is independent of Fτ and has the same law as X and hence is a Lévy
process. For instance, the first entrance time (first hitting time) of a given subset B ⊆ R is F-stopping time.
8 1.3. PRELIMINARIES ON LÉVY FLUCTUATION THEORY One elementary but useful concept, applying to all Lévy processes, is duality; it is a direct consequence of the stationary independent increments. In fact, duality can informally de-scribed as a kind of symmetry under time reversial. When a path of a Lévy process is re-versed in time, over a finite time horizon, the new path is distributionally equivalent. More precisely, for each t > 0 the processes{X(t−s)−− Xt : 0≤ s ≤ t} and {−Xs : 0≤ s ≤ t}
are equivalent, and have the same law.
An interesting direct consequence of this duality property concerns a relationship between the running supremum and the running infimum, which are defined by
¯
Xt:= sup
0≤s≤t
Xs, Xt:= inf
0≤s≤tXs.
The processes{ ¯Xt: t≥ 0} and {Xt: t≥ 0} are the key objects studied in fluctuation theory,
and will play an important role in the sequel of this thesis. We arrive at the following useful lemma.
Lemma 1.3.2. For each fixed t > 0, the pairs ( ¯Xt, ¯Xt− Xt)and (Xt− Xt,−Xt)have the same
distribution underP.
We now consider the running maximum and minimum up to τ (q), which represents an exponentially distributed time, with arbitrary parameter q > 0 (i.e., mean 1/q). It can be seen that if the Lévy process is not a compound Poisson process, then its maximums are obtained at unique times; we define ¯Gt := sup{s < t : ¯Xs = Xs} and Gt := sup{s < t :
Xs= Xs}. We are now ready to state the Wiener-Hopf decomposition (where it is noted that
the compound Poisson case should be treated slightly differently; see [68]).
Theorem 1.3.3. The Wiener-Hopf factorization. Suppose that X is any Lévy process and let τ (q)be an independent exponentially distributed random variable with parameter q > 0. Then the following statements hold.
1. The pairs
( ¯Gτ (q), ¯Xτ (q))and (τ (q)− ¯Gτ (q), ¯Xτ (q)− Xτ (q))
are independent and infinitely divisible. For any θ, ϑ ∈ R the following factorization applies: q q− iϑ + ξ(θ) =E eiϑ ¯Gτ(q)+iθ ¯Xτ(q) E eiϑGτ(q)+iθXτ(q) where the pairsE(eiϑ ¯Gτ(q)+iθ ¯Xτ(q))
andE(eiϑGτ(q)+iθXτ(q))
are the Wiener-Hopf factors. 2. When setting ϑ = 0, the Wiener-Hopf factors may be identified in terms of the Laplace
CHAPTER 1. INTRODUCTION 9 exponent κ(α, q) and ¯κ(α, q), which are defined by (for some k0> 0, and α≥ 0)
κ(α, q) :=Ee−α ¯Xτ(q) = k0exp − ∞ 0 (0,∞) 1 t e−qt− e−qt−αxP(Xt∈ dx)dt (1.1)
for the running maximum, and (for some ¯k0> 0, and α≤ 0)
¯ κ(α, q) :=Ee−αXτ(q) = ¯k0exp − ∞ 0 (−∞,0) 1 t e−qt− e−qt−αxP(Xt∈ dx)dt (1.2) for the running minimum. In addition,
κ(α, q)¯κ(−α, q) = q
q− log Ee−αX1 =:K (α, q).
Note that there are some constants in the expressions of the Wiener-Hopf factorization which are not identified (k0and ¯k0, that is). They depend on the normalization which is chosen in
the definition of local time; for more background on this issue we refer to [68].
We do not provide a proof of the above result. Instead, in order to convey the main ideas behind it, we include a rough sketch of a possible proof in a discrete-time framework (which corresponds to a random walk). We have chosen to do so, since in the continuous-time framework there are number of rather technical steps that have to be dealt with; at the same time, we believe that the main ideas behind the discrete-time counterpart provide useful insights [34].
To this end, consider the random walk Sn:=
n
i=1Yi, with the Yibeing i.i.d., distributed as
a generic random variable Y. Let ¯Snthe running maximum process
¯
Sn:= sup i∈{1,...,n}
Si;
Gn denotes the (first) epoch at which that running maximum is attained. Let T be an
(in-dependent) geometric random variable, i.e.,P(T = k) = p(1 − p)k
, for some p∈ (0, 1), and k∈ {0, 1, · · · }.
Realize that the number of maximums which are attained before time T (number of excur-sions) is a geometric random variable; it is denoted by N . It follows that both ¯ST and GT
can be written as the sum of N i.i.d. non-negative random variables. It can be showed that a geometric sum of i.i.d. random variables is infinitely divisible. Based on the above, we can
10 1.3. PRELIMINARIES ON LÉVY FLUCTUATION THEORY conclude that ¯ST and GT are infinitely divisible as well.
Furthermore, in line with the duality property of Lévy processes, we have that (ST− ¯ST, T−
GT)is independent of ( ¯ST, GT). This can be intuitively understood, as follows. First
real-ize that the geometric distribution is memoryless. Suppose now that we are told that the maximum (before time T ) is attained at a specific epoch (GT), the value of this specific issue
does not have any impact on the amount by which the process goes down between GT and
T. A similar property holds for the residual time T− GT until ‘the geometric clock expires’.
Also, it is observed that, from the duality property, ST− ¯ST has the same distribution as the
running minimum process.
After these first observations, we now include a bit of elementary algebra, leading to the identification of the Wiener-Hopf factors. To this end, we first notice that it can be verified that, with s∈ (0, 1] and α ∈ R,
E sT
eαiST = p
1− (1 − p)s EeαiY .
On the other hand, this quantity can be alternatively written as
exp − ∞ −∞ ∞ n=1 1 n(1− s n eαix)(1− p)nP(Sn ∈ dx) = exp −∞ n=1 1 n(1− s nEeαiSn )(1− p)n = exp −∞ n=1 1 n (1− p)n− (1− p)sEeαiξn = exp log p− log 1− s(1 − p) EeαiY.
Recall that we found that (ST, T )can be written as the sum of two independent terms, viz.
( ¯ST, GT)and (ST − ¯ST, T − GT)which are both infinitely divisible. As a result, it follows
that E sGT eαi ¯ST = exp − ∞ 0 ∞ n=1 1 n(1− s n eαix)(1− p)nP(Sn∈ dx) , and E sT−GT eαi(ST− ¯ST) = exp − 0 −∞ ∞ n=1 1 n(1− s n eαix)(1− p)nP(Sn∈ dx) .
substan-CHAPTER 1. INTRODUCTION 11 tial amount of technicalities. If we would ‘extrapolate’ our discrete-time findings, we would obtain the following. Let T be exponentially distributed random time with mean 1/θ, inde-pendent of the Lévy process (Xt)t. For β≥ 0 and α ∈ R, we can easily show that
Ee−βT +αiXT = ϑ
ϑ + β− log EeαiX1. (1.3)
Using the Frullani integral identity [68], we also have exp − ∞ 0 ∞ −∞ 1 t e−ϑt− e−(ϑ+β)teαix P(Xt∈ dx)dt = exp − ∞ 0 1 t e−ϑt− e−(ϑ+β)tEeαiXt dt = exp − ∞ 0 1 t e−ϑt− e−(ϑ+β−log EeαiX1)t dt = ϑ ϑ + β− log EeαiX1.
Mimicking the discrete-time setup, we obtain the same results in continuous-time frame-work: Ee−βGT+αi ¯XT = κ(ϑ + β,−αi) κ(ϑ, 0) = exp − ∞ 0 ∞ 0 1 t e−ϑt− e−(ϑ+β)teαix P(Xt∈ dx)dt , and Ee−β(T −GT)+αi(XT− ¯XT) = ¯κ(ϑ + β, αi) ¯ κ(ϑ, 0) = exp − ∞ 0 0 −∞ 1 t e−ϑt− e−(ϑ+β)teαix P(Xt∈ dx)dt . where the functions κ and ¯κare defined in Equations (1.1) and (1.2).
1.3.2
Second factorization identity
Much of our analysis relies on the Wiener-Hopf factorization and its ramifications. A second result that we use in this thesis is usually referred to as the second factorization identity and it holds for any Lévy process. It can be found in e.g. [68, p.176].
12 1.3. PRELIMINARIES ON LÉVY FLUCTUATION THEORY Let us define the first hitting time (or first passage time) as
σ(x) := inf{t ≥ 0 : Xt> x},
and let, as before, T be an exponentially distributed random variable with mean q−1. Then the following result holds for any Lévy process (Xt)t. As the proof is straightforward and
insightful, we decided to include it.
Lemma 1.3.4. For q, ¯q≥ 0, β > 0, ∞ 0 e−βxE e−qσ(x)−¯q(x−Xσ(x))1 {σ(x)<∞} dx = 1 β− ¯q 1−Ee−β ¯ XT Ee−¯q ¯XT = 1 β− ¯q 1−κ(β, q) κ(¯q, q)
where κ(α, q) was defined in Equation (1.1).
Proof. We follow the proof of [68, Exercise 6.7]. Xtis a Lévy process and hence it is
Marko-vian; in addition T is exponentially distributed. Due to the memoryless property of the exponential distribution, we have
E e−¯q( ¯XT−Xσ(x)) =E e−¯q ¯Xt and therefore E e−¯q ¯XT1 { ¯XT>x} =E e−¯q ¯XT1 {σ(x)≤T } =E e−¯qXσ(x)1 {σ(x)≤T }E e−¯q ¯XT . In addition, the first factor of the previous equation can be expressed explicitly in terms of distribution of σ(x) and T . It follows that
E e−¯qXσ(x)1 {σ(x)<T } = ∞ 0 e−qs ∞ 0 qe−q(t−s)E 1{s<t}e−¯qXsdtP(σ(x) ∈ ds) = E e−qσ(x)−¯qXσ(x)1 {σ(x)<∞} . Combining the above results gives
∞ 0 (β− ¯q)e−(β−¯q)xE e−¯q ¯XT1 { ¯XT>x} dx = ∞ 0 ∞ x (β− ¯q)e−(β−¯q)xe−¯qudP( ¯XT ∈ du)dx.
CHAPTER 1. INTRODUCTION 13 By interchanging the order of integrations we have
∞ 0 (β− ¯q)e−(β−¯q)xE e−¯q ¯XT1 { ¯XT>x} dx =E e−¯q ¯XT − E e−β ¯XT .
Theorem 1.3.3 proves the claim. 2
1.3.3
Some remarks on Wiener-Hopf factorization
The Wiener-Hopf factorization theorem provides an elegant decomposition in terms of the characteristic functions of the running maximum and the running minimum associated with the underlying Lévy process. The result, however, does not say how one should calculate each characteristic functions; preferably one would express the Wiener-Hopf factors in terms of the model primitives, i.e., the characteristic exponent ξ(·).
In fact such a decomposition is possible for specific classes of Lévy processes only. This is for instance the case if the driving Lévy process belongs to the class of spectrally one-sided Lévy processes, i.e., Xthas only negative jumps (Xt∈ S−) or only positive jumps (Xt∈ S+): in
both cases κ(α, q) can be expressed in closed-form in terms of the characteristic exponent. Let Xt∈ S−, and let Φ(β) := logEeβX1be the so-called Laplace exponent, where we define its
right inverse by Ψ(.). Then, as it turns out,
κ(α, q) = Ψ(q) Ψ(q) + α.
In case of spectrally positive case Xt ∈ S+ the decomposition is usually referred as the
(generalized version of the) Pollaczek-Khinchine formula. Then we have κ(α, q) = q
ψ(q)
ψ(q)− α q− φ(α)
where the Laplace exponent is defined by φ(α) := logEe−αX1, and ψ(.) is the inverse of φ(.).
The function ¯κ(α, q)follows from κ(α, q)¯κ(−α, q) = K (α, q) [34].
The following case can be dealt with (semi-)explicitly as well. If the jumps in one direction (either downward or upward) have a phase type distribution (which we further comment on below), whereas the jumps in the other direction are allowed to have a general distribution, the Wiener-Hopf decomposition can be performed in terms of the roots of the equation q = ξ(s); see e.g. [72, 71].
14 1.3. PRELIMINARIES ON LÉVY FLUCTUATION THEORY more explicit terms, is the class of processes which have a meromorphic Lévy exponent in the complex plane, e.g. expressed in terms of beta and digamma functions. Also for this class of Lévy processes the Wiener-Hopf factorization can be evaluated in terms of the roots of equation q = ξ(s), but now this equation has infinitely many roots [66, 67].
Any distribution can be arbitrarily well approximated by a phase type distribution [10, Thm. III.4.2]; the class of phase-type distributions is dense (in the sense of weak convergence) in the set of all probability distributions on (0,∞). If we are in a situation in which the jumps in both directions are general, we could replace those in one directions by their phase-type counterpart, leading to a Lévy process that is covered by [72, 71]. Several methods have been developed to deal with approximating a distribution on (0,∞) by a phase-type distribution, see for example [40, 57]. The approach which we used in this thesis is based on the expectation-maximization algorithm.
So far we did not discuss the impact of the ‘small jumps’ in case the driving Lévy process has infinite activity. In order to enter the setup of [72, 71], with phase type jumps in one directions and general jumps in the other direction, it is implicitly assumed that the Lévy process is of finite activity, i.e.−∞∞ Π(dx) <∞ (as these small jumps cannot be described by a compound Poisson stream of phase-type distributed jumps). This issue can be remedied as follows; focus on the situation that we wish to replace the positive jumps by a phase-type counterpart. The jump distribution on ( ,∞), for some > 0, can be approximated by a phase type distribution for any specific > 0. A Brownian motion with drift (with appropriately chosen parameters) can then compensate the jumps smaller than [10]. By picking the value of suitably small, the approximation turns out to be highly accurate. More concretely, assuming the upward jumps are approximated by the phase-type distribution Pph(dx), the parameters of the Brownian motion are calculated by the following formulas:
μ := 0 x (Π− Pph) (dx), σ2 := 0 x2(Π− Pph) (dx).
The approximation of small jumps with a deterministic drift process and Brownian motion are also frequently used in Monte Carlo simulation [10].
From the above, we conclude that if one manages to (sufficiently accurately) approximate the Lévy measure with a phase type distribution, the factors in the Wiener-Hopf decompo-sition can be evaluated. As a consequence, we have the Laplace transforms of the running maximum and minimum. The next step is to invert these, in order to evaluate the
distri-CHAPTER 1. INTRODUCTION 15 butions of the running maximum and minimum. There are several methods developed for performing Laplace (and/or Fourier) inversion. Most of the Laplace inversion techniques are based on the well-known Poisson summation formula (PSF) [1, 38]. ThePSFis given by the following formula, for any v∈ [0, 1) and any ‘damping factor’ a ∈ R:
∞ k=−∞ ˆ f (a + 2πi(k + v)) = ∞ k=0 e−ake−2πikvf (k)
where ˆf is the Laplace (Fourier) transform of the function f (x). The right hand side of this equation is a discrete Fourier transform which can be computed efficiently by the well-known fast Fourier transform algorithm, obviously provided that one can evaluate the left-hand side of the equation. The method which is developed by den Iseger [36] approxi-mates the infinite summation with a finite sum at appropriately chosen points and weight. This technique has been extensively tested, and has shown to be able to calculate Laplace and Fourier inversion transform fast and accurately. The technique can be adapted to per-form the inversion of non-smooth functions and even functions with singularities. It is also remarked that the extension of the method to multi-dimensional mixed Laplace/Fourier in-version is straightforward; for details of the implementation, as well as a series of extensions, we refer to [36].
1.4
Preliminaries on Markov fluid models
The topics of fluctuation theory and queueing theory are intimately related; e.g. the steady-state distribution of the workload in a queueing system can often be translated in terms of the probability of an associated ‘free’ (i.e., not truncated at 0) process attaining a given set. In this sense, techniques developed in the context of fluctuation theory for Lévy processes, can be used in the context of Lévy-driven queues as well [34, 83]. In this section we leave the setting of Lévy processes though, focusing on queues with Markov fluid input. We provide the preliminaries necessary for the last chapter of this thesis.
Queueing theory studies the evolution of a storage process, which can be in terms of either (discrete) customers or workload. A queue is characterized by an arrival process, the distri-bution of the service requirements, and a service discipline. In general a queue can have one or more servers (where servers can be interpreted as internet servers, cashiers in shops, etc.), which process the clients’ jobs. Often the arrival process and service times are uncertain, and therefore random objects are used to model these.
16 1.4. PRELIMINARIES ON MARKOV FLUID MODELS There is a vast body of literature on the mathematical modeling of all sorts of queues. A large subclass of these models can be summarized by the notation A/B/n/m− S, which characterization is due to Kendall in the 1950s. Here A and B correspond to, respectively, the distributional properties of the arrival process and service requirements. The number of servers and the maximum number of jobs which can wait in the system (i.e., in its buffer) are indicated by n and m respectively. Finally, the S specifies the service discipline; it can be for instance first-come-first serve or processor-sharing.
The model we consider in this thesis, does not fit in the Kendall notation. We consider a fluid model, in which a reservoir is fed by a continuous traffic stream [65, 83, 4, 63, 30]. The server may be considered as the output flow of the reservoir; the output rate is usually considered constant but may also be stochastic [74, 75, 89, 42, 20]. We give a more detailed description below.
1.4.1
Markov fluid model
Consider the following fluid reservoir (or buffer), where the amount of fluid in the reservoir at time t is denoted by Ct. Let (Xt)tdenote the so-called background process, which is
as-sumed to be an irreducible continuous Markov process; this background process models the stochasticity of the input flow into the reservoir. We assume that Xthas a finite state space
N ⊂ N, i.e., Xtattains values in the setN = {1, 2, · · · , N}.
Now the content of the reservoir is driven by (Xt)t, in the sense that the input rate into the
reservoir is ri when the process (Xt)tis in state i ∈ N , unless the reservoir is empty and
the net input flow rate is negative (as in that situation the reservoir remains empty). The content of the reservoir is evidently stochastic, and its dynamics are given by the following differential equation: dCt dt = ⎧ ⎪ ⎨ ⎪ ⎩ max(ri, 0) if Ct= 0, ri if Ct> 0. (1.4) The above formula tacitly assumed that the buffer has infinite capacity; as an aside, we note that if there is a finite buffer B > 0, the following equation holds:
dCt dt = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ max(ri, 0) if Ct= 0, min(ri, 0) if Ct= B, ri if Ct∈ (0, B).
CHAPTER 1. INTRODUCTION 17 It is assumed that at least one of the net flow rates is positive, as otherwise the model is trivial (in that Ct= 0all the time). However, when the buffer capacity is infinite, the requirement
that the queue be stable means that we have to impose the condition
i∈N
piri< 0,
where piis the stationary probability of Xtbeing in state i∈ N ; if this condition would not
be met, the queue will grow to infinity [65].
We define by Fi(y, t)the probability of the content is at most y and the background process
Xtis in state i. In other words,
Fi(y, t) :=P[Xt= i, Ct≤ y], i∈ N , y ≥ 0.
In addition, the background process is a continuous-time Markov chain with generator ma-trixQ = [qij]such that
P[X(t + h) = j|X(t) = i] = qijh + o(h),
P[X(t + h) = i|X(t) = i] = 1 + qiih + o(h). (1.5)
where qij ≥ 0 if i = j and qii =−
j=iqij,for i∈ N .
Consider Fi(y, t + h), i.e., the probability of the background process being in state i and
the buffer content is being at most y at time t + h, where h ‘is infinitesimally small’. Using Equations (1.5), Fi(y, t + h) can be expressed in terms of Fj(y, t). Let y > 0, h > 0 and
y− rih > 0for all i∈ N . Then,
Fi(y, t + h) = (1 + qiih)Fi(y− rih, t) + h
j=i
qjiFj(y− rjh, t) + o(h). (1.6)
It is elementary to rewrite the above equation to
(Fi(y, t + h)− Fi(y)) + (Fi(y)− Fi(y− rih, t)) = h
j=i
qjiFj(y− rjh, t) + o(h).
18 1.4. PRELIMINARIES ON MARKOV FLUID MODELS equation: ∂Fi(y, t) ∂t = j∈N qjiFi(y, t)− ∂Fi(y, t) ∂y ri.
This equation (which can be regarded as a Kolmogorov forward equation) [65, 60, 77] can be written in compact matrix form; then it becomes
F (y, t)
∂t =F (y, t)Q −
F (y, t)
∂y R, y > 0, (1.7) whereF (y, t) = (F1(y, t),· · · , FN(y, t))T andR = diag(r1,· · · , rN). When the stability
con-dition is satisfied, then it can be shown that there exists a stationary distribution of (Xt, Ct);
i.e., the partial derivative with respect to time vanishes, and we have that, as t→ ∞,
F (y)Q = F (y)
∂y R (1.8)
whereF (y) denotes the corresponding stationary distribution of the fluid level in the reser-voir.
Evidently, (1.8) alone does not uniquely specify the stationary distribution; a set of coefficient still needs to be specified. More specifically, the solution of (1.8) has the form
F (y) =
N
j=1
cjeξjyV(j), (1.9)
where (ξj,V(j))are the eigenvalues and corresponding eigenvectors of the matrixR−1QT,
and cj are constants which have to be determined by imposing boundary conditions [69].
Particularly, when the buffer size is infinitely large, if Re(ξj) > 0, the corresponding
coeffi-cient cj has to be zero, as otherwise the probability F (y) cannot be bounded between 0 and
1. We divide the states into two separate sets, in terms of their net input rates; we define N+≡ {i ∈ N | r
i > 0} and N−≡ {i ∈ N | ri< 0} (assuming for ease that there is not i such
that ri= 0). As a consequence, we have the following boundary condition corresponding to
an empty reservoir:
Fi(0) = 0, i∈ N+;
CHAPTER 1. INTRODUCTION 19 buffer is finite then the following boundary conditions have to hold:
Fi(B−) = pi, for i∈ N−
where B is the buffer size. It turns out that these condition yield as many conditions as there are unknowns cj (viz. N ). An important role in this argument is played by the property
that, if the stability condition is satisfied, the number of eigenvalues with negative real part equals|N+|, there is one zero eigenvalue and the other eigenvalues have positive real part.
1.4.2
Workload-dependent service rate models
In the setting above, the net input rates and the transition matrix did not depend on the current buffer content (i.e., the fluid amount in the reservoir). We now consider situations in which we depart from this framework.
A first scenario is that in which the rates riand the transition matrixQ are functions of the
current buffer content. This class of models is, for obvious reasons, sometimes referred to as feedback models. We primarily consider the case in which there are (finitely many) levels such that between two subsequent levels ri andQ are constant [89, 42]. We consider a second
scenario as well, in which the rates ridepend on the direction in which the level is crossed;
we call this scenario an hysteresis-type model [74, 75].
First we consider models in which ri and Q are locally constant between specified levels,
and do not depend on the direction in which the level is crossed. Suppose there are K + 1 buffer levels such that
0 = B0≤ B1≤ · · · ≤ BK−1≤ BK ≤ ∞
and, as a consequence, K buffer regimes; when y ∈ (BK−1, BK)the flow rates are given
byRk = diag(rk1,· · · , rNk)where k ∈ {1, 2, · · · , K}. We remark that this setup can be even
extended to a continuous counterpart, in which riandQ depend on the current buffer level
in a continuous fashion [89]; here we only consider the case that flow rates are piecewise constant functions of buffer content.
Assume that allRk matrices are invertible. As a consequence, the steady-state probability distribution of buffer content in each regime k is Fk
i(y)follows from a Kolmogorov equation
in the spirit of Equation (1.8) (withR and Q replaced by their ‘local counterparts’), and the (local) solutions are given by an expression in the spirit of Equation (1.9). To complete the
20 1.5. OUTLINE OF THESIS, CONTRIBUTIONS solution we need to align the solutions corresponding to the individual regimes. This is done by imposing the appropriate boundary conditions at each buffer threshold. It requires a substantial amount of administration to identify all boundary conditions, and to verify that their number equals the number of unknown coefficients.
We now focus on the hysteresis-type model. For ease, we now assume that the buffer is infi-nite and there are only two thresholds, the lower level and upper level which are indicated by B and Bu, respectively (where obviously B ≤ Bu). In addition, for ease we assume that
in each background process state, the net input rate rate can take two values (i.e., a higher and the lower value). Introduce an indicator process I(·) ∈ {+, −}, corresponding to the two regimes, i.e., those in which the higher and lower net input rates apply.
The process is assumed to evolve as follows. It starts with an empty buffer and the indicator is +. The indicator stays + as long as the buffer content remains below Bu. At the moment
that buffer content reaches the level Bu(and the background process is in state i), the
indi-cator changes from into− (while the background process remains in i). Then the indicator remains− as long as the buffer content is above the level B ; then it changes to− again
(where the background process does not change). The process continues in this fashion. With techniques similar to those explained above, we can evaluate the steady-state proba-bilities
Fi−(x) :=P(I = −, X = i, C ≤ x), Fi+(x) :=P(I = +, X = i, C ≤ x).
Finding these is again a matter of setting up differential equations, and imposing the appro-priate boundary conditions. For a detailed analysis we refer to [74].
1.5
Outline of thesis, contributions
The primary objective of this thesis is to contribute to the development of computational techniques in queueing and fluctuation theory. Essentially, three types of techniques will be explored.
• In the first place, we systematically validate a (one- and two-dimensional) Laplace and Fourier inversion algorithm. Our approach is based on an algorithm proposed by den Iseger [36], and several variants thereof, which essentially rely on the Poisson summa-tion formula. We do so in the context of Lévy fluctuasumma-tion theory, aiming at numerically
CHAPTER 1. INTRODUCTION 21 evaluating the probability distribution of the running maximum process. While this approach has a variety of potential applications, we primarily focus on applying it to price specific exotic options, viz. the so-called lookback option.
• The second technique we present in this thesis is importance sampling. This technique aims at reducing the variance of simulation-based estimators of performance mea-sures. For instance, when estimating rare event probabilities, straightforward simu-lation is extremely time consuming, as many paths need to be generated in order to obtain an estimate with low variance. The idea behind importance sampling is to gen-erate paths under an alternative measure, which is chosen such that the event under study is not rare anymore. Correcting the simulation output by a likelihood ratio, the resulting estimator is unbiased. The challenge is to find a good new measure, for which the variance of the estimator (provably) reduces.
• In the third place, considering a node in a communication network, we assess the tradeoff between the quality of service and the energy consumption. Different service scenarios are considered; in each scenario the service speed is determined in a specific way by the evolution of the buffer occupancy. We use techniques from optimization to minimize a cost function (encompassing quality of service and energy consumption), where the parameter space covers all feasible service strategies. The optimization rou-tine is based on simulated annealing in combination with classical Newton-Raphson-type algorithms, in the sense that we identify by simulated annealing the initial point of the Newton-Raphson optimization routine (which locally finds the optimum). In more detail, the contributions of the individual chapters can be summarized as follows. There are four studies, the first three focusing on fluctuation-theoretic aspects in a Lévy-driven system, and the last evaluating the tradeoff between quality of service and energy consumption in a queue with fluid input. We now detail the specific contributions.
Chapter 2 presents a framework for numerical computations in fluctuation theory for Lévy processes. More specifically, with as before ¯Xt := sup0≤s≤tXsdenoting the running
max-imum of the Lévy process Xt, the aim is to evaluateP( ¯Xt ≤ x) for t, x > 0. We do so
by approximating the Lévy process under consideration by another Lévy process for which the double transformEe−α ¯Xτ(q)
is known, with τ (q) an exponentially distributed random variable with mean 1/q; then we use a fast and highly accurate Laplace inversion technique (of almost machine precision) to obtain the distribution of ¯Xt. A broad range of examples
22 1.5. OUTLINE OF THESIS, CONTRIBUTIONS In Chapter 3 our objective is to compute the prices (and corresponding sensitivities, known as Greeks) of lookback options driven by Lévy processes. In this setup, the risk neutral evolu-tion of the stock price, say St, is given by S0eXt, with S0the initial price and Xtrepresenting
a Lévy process. Lookback options prices are functions of the stock price ST at the maturity
time T and the running maximum ¯ST := sup0≤t≤TSt, and as a consequence the
Wiener-Hopf decomposition provides us with all probabilistic information needed to evaluate these prices. To overcome the complication that in general only an implicit form of the Wiener-Hopf factors is available, we follow the same approach as in Chapter 2: we approximate the Lévy process under consideration by an appropriately chosen other Lévy process for which the double transform Ee−α ¯Xτ(q) is known; as before, τ (q) is an exponentially distributed
random variable with mean 1/q. The second step is to write the transform of the lookback option prices in terms of this double transform. Finally, we use state-of-the-art numerical inversion techniques to compute the prices and Greeks (i.e., sensitivities with respect to ini-tial price S0and maturity time T ); these rely on the techniques featuring in Chapter 2. We
test our procedure for a broad range of relevant Lévy processes, including a number of ‘tra-ditional’ models (Black-Scholes, Merton) and more recently proposed models (CGMY and Beta processes), showing excellent performance in terms of speed and accuracy. This has been submitted for publication to Journal of Computational Finance [7].
In Chapter 4 the focus is on numerical techniques to evaluate rare-event probabilities in a Lévy setting. We analyze the tail asymptotics corresponding to the all-time maximum value attained by a Lévy process with negative drift. This chapter has two main contributions: a short and elementary proof of these asymptotics, and an importance sampling algorithm to estimate the rare-event probabilities under consideration. This chapter is based on [6] which has been accepted for publication in Statistics and Probability Letters.
Finally, in Chapter 5 considers the tradeoff between quality of service and capacity cost in communication networks. More specifically, we develop techniques for analyzing and op-timizing energy management in multi-core servers with so-called speed scaling capabilities (i.e., the service speed can be adjusted based on the current buffer occupancy, or the evolu-tion of the buffer occupancy in the recent past). Our framework incorporates the processor’s dynamic power, but it also accounts for other intricate and relevant power features such as the static (leakage) power and switching overhead between speed levels. Using stochas-tic fluid models to capture traffic burst dynamics, the chapter proposes and studies differ-ent strategies for adapting the multi-core processor speeds based on the observable buffer content, so as to optimize objective functions that balance energy consumption and
perfor-CHAPTER 1. INTRODUCTION 23 mance. The strategies can be non-hysteretic (i.e., the processor speed depends on current buffer level relative to the buffer thresholds) or hysteretic (i.e., it matters in which direction the buffer thresholds are crossed). It is shown that, under rather general conditions, strate-gies which use more threshold levels are more efficient with respect to power consumption; however, most of the efficiency gain is achieved with 1 or 2 thresholds only. In addition, the optimal power consumptions of the different strategies are only very mildly sensitive to perturbations in the input parameters, implying the highly advantageous property that the performance is robust to estimation errors in the system’s input traffic parameters. This chapter has appeared as [9], and a short version as [8].
Our objective is that all chapters are self-contained, i.e., they can be read separately. As a consequence, there will be some inevitable amount of overlap between these chapters. We also remark that we have pursued a maximum level of uniformity regarding the notation used throughout the thesis, and that this is also in line with the notation introduced in the present chapter.
Chapter
2
Numerical techniques in Lévy
fluctuation theory
In this chapter we discuss a numerical technique to fast and accurately evaluate the distri-bution of the running supremum as attained by a Lévy process.
2.1
Introduction
As explained in Chapter 1, owing to their wide applicability and their attractive mathemati-cal properties, Lévy processes play an important role in applied probability. In mathematimathemati-cal terms, they are characterized as processes with stationary and independent increments, and, as such, the class of Lévy processes covers e.g. Brownian motion and (compound) Poisson processes (but is substantially broader; for instance processes with infinitely many jumps in finite time intervals belong to this class as well). Over the past decades Lévy processes have found widespread use in various application domains. More specifically, they are in-tensively studied in both mathematical finance and operations research, see, among many other sources, for instance [10, 31, 35].
With Xt denoting the Lévy process (assuming X0 = 0), a substantial research effort
con-centrates on analyzing probabilistic properties of the so-called running maximum process ¯
Xt:= sup0≤s≤tXs.More particularly, one wishes to determine the probabilityP( ¯Xt≤ x) for
t, x > 0, or alternatively the corresponding density. The branch of research focusing on this type of problems is commonly known as fluctuation theory [21, 68, 83].
26 2.1. INTRODUCTION As mentioned in Chapter 1, a Lévy process is characterized by its Lévy exponent logEeisX1
, which is a necessarily of the form
logEeisX1 = isd−1
2s 2 σ2+ ∞ −∞(e isx− 1 − isx1 {|x|<1})Π(dx), (2.1)
where d∈ R, σ ≥ 0, and the spectral measure Π(·), concentrated on R \ {0}, satisfies
R
min{x2, 1}Π(dx) < ∞. The triplet (d, σ2
, Π)is usually referred to as the characteristic triplet, as it uniquely defines the Lévy process [21, Ch. I, Thm. 1]. The three terms in the right-hand side of the representation (2.1) are, for obvious reasons, often called the (deterministic) drift term, the Brownian term, and the jump term. Special cases of Lévy processes are deterministic drifts (only a drift term) and Brownian motions (only a Brownian term). The class of Lévy processes also contains compound Poisson processes; then we just have the jump-term (and the first term as well in case a deterministic drift is present as well), and in addition there should be a well-defined arrival rate (which requires that −∞∞ Π(dx) < ∞). The class is wider though, as it also includes processes with infinitely many jumps in a finite amount of time (usually referred to as ‘small jumps’); this happens in case−∞∞ Π(dx) =∞.
In principle the distribution of ¯Xt is fully specified through the so-called Wiener-Hopf
de-composition, see e.g. [68, Ch. 6]. It implies that, with τ (q) denoting an exponential random variable with mean 1/q that is independent of the Lévy process Xt,
κ(α, q) :=Ee−α ¯Xτ(q)= k 0exp − ∞ 0 (0,∞) 1 t e−t− e−qt−αxP(Xt∈ dx)dt , (2.2)
where k0 is a normalizing constant. From a practical standpoint, the use of this
characteri-zation is limited, as it provides us with the double transform ofP( ¯Xt∈ dx) — realize that
1 q · Ee −α ¯Xτ(q) = ∞ 0 e−qt ∞ 0 e−αxP( ¯Xt∈ dx)dt, (2.3)
which in general cannot be inverted explicitly.
The above entails that, in order to get numerical values for the densityP( ¯Xt ∈ dx) or the
distribution functionP( ¯Xt ≤ x), one option is to (i) first evaluate the double integral (2.2)
