On the power of formal contract description for quantitative purposes: an example with control variates

Vivien Begot 2012-06-15
As my first subject (and the first subject on LexiFi's quantitative blog), I would like to focus on something I particularly appreciate at LexiFi: our formal contract description, that allows one to treat quantitative topics in a surprisingly elegant and easily generalisable way.
I will focus on a particular topic: control variates methods. In the first section, I will introduce the underlying theory. Then, I will explain its use in mathematical finance. To finish, I will try to show the interest of a formal contract description for this particular topic, writing a generic control variable selection algorithm.

Theoretical introduction

Monte Carlo and variance reduction

Monte Carlo method is based on the approximation of an expectation using the law of large numbers. Let $X$ be a random variable, $(X_i)_{i\in I}$ a family of i.i.d random variables with the same law as $X$. Let $\bar X_n = \frac{1}{n}\sum_{i=1}^n X_i$ $\bar X_n$ converges in probability to $\mathbb{E}[X]$, and $\mathbb{E}[\mathbb{E}[X] - \bar X_n ] = 0$ (i.e. $\bar X_n$ is an unbiased estimator of $\mathbb{E}[X]$). Moreover, using the central limit theorem, $\sqrt n \frac{\bar X_n - \mathbb{E}[X]}{Var(X)}$ converges in law to a centered gaussian law with variance one. It shows that the variance controls the convergence rate of our estimator. The goal of variance reduction is then to find a variable $\tilde X$ such that $\mathbb{E}[\tilde X] = \mathbb{E}[X]$ and $Var(\tilde X) < Var(X)$. Control variates is a particular method to achieve this goal.

Control variates method

The idea of this method is to take a linear combination of others random variables (with null expected value), and to add it to the initial one. The new variable will have the same expected value as the first one, but there is no reason for them to have the same variance.

Let $Y = (Y^1, ..., Y^d)^T$ be a $d$-dimensional real valued random variable, with mean $\mu_Y$. We define: $\tilde X_\beta = X - \beta^T(Y - \mu_y)$ Where $\beta = (\beta ^ 1, ..., \beta ^ d) ^ T$. Then, \begin{aligned} Var(\tilde X_\beta) & = Var(X) + \beta ^ T \Sigma_Y \beta - 2 \beta ^ T \Sigma_{XY}\\ \mathbb{E} [ \tilde X_\beta ] & = {\mathbb{E}} [ X ] \end{aligned} Where $\Sigma_Y$ is is the covariance matrix of $Y$, and $(\Sigma^i_{XY})_{i \in [1; d]} = \left(Cov(X, Y^i)\right)_{i \in [1; d]}$. For more convenience, we assume that $\Sigma_Y$ is positive-definite (this is not too restrictive as it is at least positive semi-definite).
At this stage, we want to find $\beta^*$ such that $\tilde X := \tilde X_{\beta ^ *}$ has the lowest possible variance, i.e. $Var(\tilde X) = min\{ Var(\tilde X_\beta); \beta \in \mathbb R ^ d \}$ Such a $\beta$ exists because the variance is a positive second order polynomial of $\beta$. Let $v (\beta) = Var (X_\beta)$. $\beta^*$ verifies:
\begin{aligned} \nabla v (\beta ^ *) & = 2\Sigma_Y \beta ^ * - 2\Sigma_{XY} = 0\\ \Rightarrow \beta ^ * & = \Sigma_Y^{-1}\Sigma_{XY} \end{aligned} where $\nabla$ is the gradient operator. Moreover, $Var (\tilde X) = Var(X) - \Sigma_{XY} ^ T \Sigma_Y ^ {-1} \Sigma_{XY} \leq Var(X)$

Remark: if we don't restrict $\Sigma_Y$ to be positive-definite, it may be not invertible. However, there exists at least one $\beta$ that minimises $v$, and it can be found by replacing $\Sigma_Y^{-1}$ by the pseudo-inverse of $\Sigma_Y$. Indeed, $\Sigma_Y^+\Sigma_{XY}$ will minimize the function $\nabla v (\nabla v)^T$, and we know that the minimum is 0. Moreover, $\Sigma_Y ^ {-1} = \Sigma_Y^+$ when $\Sigma_Y$ is invertible, then it is a good idea to always use pseudo-inverse to find $\beta^*$.

That's it, we know all what we need about control variates method! But now, the real difficulty is to find the good $Y$, and this is where the formal contract description takes all its meaning. Before showing the interest of formal contract description, I suggest to make a brief section on the application to mathematical finance, so be patient, I keep the best for the end.

Application to mathematical finance

We recall that (under the classical mathematical finance assumptions) the price of a product is given by the present value of our numeraire multiplied by the expectation of the cash-flows "discounted" with the numeraire, under the risk-free probability associated to this numeraire.
Mathematically, let $S(t) = \big(S^1(t), ..., S^n(t)\big)$ be the value at t of $n$ tradeable assets. We assume that the market generated by $S$ is complete, and we work on it. Let $N$ be one of these assets, and let $V (t)$ be the value at $t$ of the product that pays at $t_i$ $\phi_i \big( \{(u, S(u)); 0 \leq u \leq t_i \}\big)$ on the schedule $(t_1, ..., t_m)$. We have $V(t) = N(t) \mathbb E^N \left[\sum_{i = 1} ^ m \phi_i \big( \{(u, S(u)); 0 \leq u \leq t_i \}\big) N(t_i)^{-1} | \mathcal F_t\right]$ where $\mathbb E^N$ is the expectation under the risk-free probability associated to $N$, and $\{\mathcal F_t; 0 \leq t\}$ the filtration generated by $\{S(u); 0 \leq u \leq t\}$ completed to satisfy the "usual conditions" (don't worry, it is not that important for this article).

This expectation can be computed by a Monte Carlo method, and we may want to apply control variates method to it. Here we have to find our $Y$. A good choice is to take simple financial products for which we have closed forms, like discounted forward, or simple european calls. We note that using notation of the previous section, we have $X = \sum_{i = 1} ^ m \phi_i \big( \{(u, S(u)); 0 \leq u \leq t_i \}\big) N(t_i)^{-1} \big| \mathcal F_t$. For more convenience for the following, we will take $Y$ as a matrix and not as a vector (it is easy to transform the matrix in vector). Thus, a simple (and not so bad) choice for $Y$ could be $Y = (Y_{ij})_{(i, j) \in [1; n] \times [1; m]}$ where $Y_{i, j} = S^i(t_j) N(t_j) ^ {-1}$

We can even improve this $Y$ by restricting ourself to contracts that just observe the underlyings on a finite number of dates, i.e. the i-th cash-flow becomes $\psi_i \left(S(u^i_1), S(u^i_2), ... ,S(t_i)\right)$. In this case, we take $Y_{i, j, k} = S^i(u^j_k) N(u^j_k) ^ {-1}$ We will see in the next section the ease of find these dates with formal contract description. But wait, I tried to cheat and "forgot" to speak about european call. Well, the call control-variables will be of the form $Y = \left(S^i(u) - K\right)^+ N(u) ^ {-1}$ and will (of course) be useful when some cash flows are close to call's one. To this effect we need to know the $\psi_i$, and if we want our algorithm to be generic, I do not see any other way than using formal contract description.

The power of formal contract description

I will start this section by briefly explaining what I mean by "formal contract description". Actually, and to make it easy, that's just the idea of creating a language to describe contracts. To be clear, I will introduce a simple grammar that will be the framework for this section:

Contract $\rightarrow$ Contract AND Cash_flow
Contract $\rightarrow$ Cash_flow

Cash_flow $\rightarrow$ PAY Amount AT Date

Amount $\rightarrow$ (Binop Amount Amount)
Amount $\rightarrow$ (Unop Amount)
Amount $\rightarrow$ (UNDERLYING String FIXED_AT Date)
Amount $\rightarrow$ (FIXED_FLOW Float)

Unop $\rightarrow$ EXP | LOG | SQRT | ...

Binop $\rightarrow$ ADD | SUB | MUL | DIV | MAX | MIN | ...

For example, let us consider a simple european call on S, with option maturity $t\_ex$, that pays at $t\_pay$, with strike $K$. This contract is written in our grammar:

PAY (MAX 0 (SUB (UNDERLYING S FIXED_AT $t\_ex$) (FIXED_FLOW $K$))) AT $t\_pay$

With such a contract representation, it is easy in a first step to get at least all interesting dates (and to build the first $Y$ on discounted forwards); it just consists in a naive tree traversal, where on each node you add the list of dates you get from the sub-trees, and add it to the date of this node if there is one.
Let's write some OCaml to have fun:

In a first step, we have to represent our grammar. The best way to do it is to use algebraic data type:

 1 2 3 4 5 6 7 8 9 10 11 12 13    type binop = Add | Sub | Mul | Div | Max | Min (* | ... *)   type unop = Exp | Log | Sqrt (* | ... *)     type amount =     | Binop of binop * amount * amount     | Unop of unop * amount     | Fixed_underlying of string * date     | Fixed_flow of float     type cash_flow = amount * date     type contract = All of contract * cash_flow | Single_flow of cash_flow   

Then, we define the call:

 1 2 3 4 5    let call k s t_ex t_pay =     let ul = Fixed_underlying ("Eurostoxx", t_ex) in     let flow = Binop (Sub, ul, Fixed_flow k) in     Single_flow (((Binop (Max, Fixed_flow 0., flow))), t_pay);;   

and we appreciate the proximity between the grammar production rules (resp. an element on this grammar) and the corresponding algebraic data type (resp. an element of this type).

The function that collects the dates will thus be:

 1 2 3 4 5 6 7 8 9 10    let rec find_dates_in_flow = function     | Binop (_, f1, f2) -> (find_dates_in_flow f1) @ (find_dates_in_flow f2)     | Unop (_, f) -> find_dates_in_flow f     | Fixed_underlying (_, t) -> [t]     | Fixed_flow _ -> []     let rec find_dates = function     | All (c, (f, t)) -> t :: ((find_dates c) @ (find_dates_in_flow f))     | Single_flow (f, t) -> t :: (find_dates_in_flow f)   

Now that we are comfortable with OCaml, why don't we think of an algorithm that detects close-to-call payoffs ? OK, let's go. A first remark is that when we are on a call, we are in a subtree of a max. A second one is that the strike is determined by something like Binop (Sub, o1, o2) where o1 (resp. o2) has a child that is a fixed underlying and the child of o2 (resp. o1) is a simple computation on fixed flows (and determines the strike). These two remarks are a good start point to build a heuristic to find close-to-call payoffs, and I propose the following function, that returns a list of (string * date * float), where string is the underlying identifier, date the maturity of the call and float the strike:

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62    let apply_unop = function     | Exp -> exp     | Log -> log     | Sqrt -> sqrt   ;;     let apply_binop = function     | Add -> (+.)     | Sub -> (-.)     | Mul -> ( *.)     | Div -> (/.)     | Max -> max     | Min -> min   ;;     let rec find_underlying = function     | Binop (_, f0, f1) -> (find_underlying f0) @ (find_underlying f1)     | Unop (_, f) -> find_underlying f     | Fixed_underlying (s, t) -> [s, t]     | Fixed_flow _ -> []   ;;     let rec find_float = function     | Binop (b, f0, f1) ->         begin match find_float f0, find_float f1 with         | Some f0, Some f1 -> Some (apply_binop b f0 f1)         | _ -> None         end     | Unop (u, f) ->         begin match find_float f with         | Some f -> Some (apply_unop u f)         | None -> None         end     | Fixed_underlying _ -> None     | Fixed_flow f -> Some f     let rec find_calls_flow in_call = function     | Binop (binop, f0, f1) ->         if binop = Max || binop = Min then           (find_calls_flow true f0) @ (find_calls_flow true f1)         else if binop = Sub && in_call then           let test = function             | Some k, (_ :: _ as l) -> List.map (fun (s, t) -> (s, t, k)) l             | _ -> []           in           let l1 = test (find_float f0, find_underlying f1) in           if l1 = [] then             let l2 = test (find_float f1, find_underlying f0) in             if l2 = [] then               (find_calls_flow in_call f0) @ (find_calls_flow in_call f1)             else l2           else l1         else (find_calls_flow in_call f0) @ (find_calls_flow in_call f1)     | Unop (_, f) -> find_calls_flow in_call f     | Fixed_underlying _ -> []     | Fixed_flow _ -> []     let rec find_calls = function     | All (c1, (f, _)) -> (find_calls c1) @ (find_calls_flow false f)     | Single_flow (f, _) -> find_calls_flow false f   ;;   
You can now try this at home, and you should have these results:

 1 2 3 4 5 6    # find_dates (call 2500. "Eurostoxx" 2013-01-01 2013-01-03);;   date list = [2013-01-03; 2013-01-01]     # find_call (call 2500. "Eurostoxx" 2013-01-01 2013-01-03);;   (string * date * float) list = [("Eurostoxx", 2013-01-01, 2500.)]   

Harder: a call on the average value of "Eurostoxx" on two dates: $max (0, 0.5 * (S_{t_1} + S_{t_2}$) - K):

 1 2 3 4 5 6 7 8    let average_call =     let ul t = Fixed_underlying ("Eurostoxx", t) in     let sum = Binop (Add, ul 2012-06-01, ul 2013-01-01) in     let average = Binop (Mul, Fixed_flow 0.5, sum) in     let flow = Binop (Sub, average, Fixed_flow 3500.) in     let payoff = Binop (Max, Fixed_flow 0., flow) in     Single_flow (payoff, 2013-01-01)   ;;
 1 2 3 4 5 6 7    # find_dates average_call;;   date list = [2013-01-01; 2012-06-01; 2013-01-01]     # find_calls average_call;;   (string * date * float) list =   [("Eurostoxx", 2012-06-01, 3500.); ("Eurostoxx", 2013-01-01, 3500.)]   

We see on this example that our simple heuristic even works for asian options. Obviously, these simple study cases can be easily enriched, what would have not been interesting in this short presentation.

Conclusion

What I tried to show with this light example is that the use of formal contract description is a powerful tool to use in a generic way quantitative methods that seemed to be not easily generalisable prima facie. Even if some part of our test can be made without our approach (I have no doubt that anybody serious can easily get the dates of his contracts), I think that it allows one to have nice algorithm for cases that could be easily treated and to treat some non-trivial cases in a generic way, while nicely handling (always generically) many other topics, like Monte Carlo and PDE pricing ocaml generation, risk analysis, fixing management, etc.
In LexiFi Apropos, you can appreciate a richest language, with finest control variable detection implemented, and other kind of control variables handled (like down-in or up-in digital call).
I hope you appreciated reading this article, and that I succeeded in showing you (even just a little) the interest of our approach for quantitative purposes.

LexiFi • 892 rue Yves Kermen • F-92100 Boulogne-Billancourt • France