two anniversary dates. PDE (2.12) is solved using the Crank-Nicolson finite differences methods. [60], [99] used the scheme for pricing GMWB with discrete optimal withdrawals. Of course, if the volatility and/or interest rates are stochastic, then one needs to add extra dimensions to the PDE.

Applying the jump condition (2.10) to the solution at $t_{N-1}^+$ we obtain the solution at $t_{N-1}^-$ from which further backward time stepping gives us solution at $t_{N-2}^+$, and so on. The numerical algorithm takes the following key steps:

1. Generate a finite grid for the account value A and benefit base G, i.e. $A_0 < A_1 < ... < A_J$ and $0 = G_0 < G_1 < ... < G_K$.

2. At $t_N$, define the final condition for each note point $(A_j, G_k)$, $j = 1, 2, ..., J$ and $k = 1, 2, ..., K$ to get $\Phi(t_n^-, A, G)$ and the boundary conditions (2.17)-(2.18) for $A_{min}$ and $A_{max}$ for each potential $G_{k \in \{1, 2, ..., K\}}$.

3. For each potential benefit $G_k$, $k = 1, 2, ..., K$, solve the PDE using the Crank-Nicolson finite differences scheme to obtain $\Phi(t_{N-1}^+, A, G)$.

4. Apply the jump condition (2.10) by performing a linear search of the withdrawal amount $\gamma_{N-1}^*$ that gives the maximum $\Phi(t_{N-1}^-, A, G)$. In general, this involves a two-dimensional interpolation in $(A, G)$ since the $h^G(A, G, \cdot)$ and $h^G(A, G, \cdot)$ do not necessarily fall in the grid nodes.

5. Repeat (3) and (4) for $t = t_{N-2}, t_{N-3}, ..., t_1$.

6. Evaluate Equation (2.12) for the backward time step $t_1$ to $t_0$ to obtain solution $\Phi(t_0, A, G)$ at $A_0$ and $G_0$.

We can add more complexity to the model, for example by incorporating stochastic interest rates or stochastic volatility. In this case, the dimension of the pricing PDE (2.12). We can also add more constant path-wise state variables that evolve only at the anniversary (tax base, extra benefit base, etc...), which will affect only the jump condition, i.e. the search of the optimum will have to be performed based on the new variables as well.

The finite difference scheme will be discussed in more detail in the next section.

## Description of the finite differences scheme

Recall that the value function $\Phi$ satisfies the following recursion

$$ \Phi(t_n^+, A, G) = \mathbb{E}_{t_n}^{\Phi} \left[ B_{n, n+1} \Phi(t_{n+1}^-, A, G, G_{n+1}) | A, G \right], \quad (2.13) $$

$$ \Phi(t_n^-, A, G) = \max_{\gamma_n \in \mathcal{A}_n} \left( p_n \bar{f}_n(A, G, \gamma_n) + p_{n-1} q_n D(t_n^+, A, G) \right. \\ \left. \qquad + \Phi(t_n^+, h^A(A, G, \gamma_n), h^G(A, G, \gamma_n)) \right). \quad (2.14) $$

Within each time interval $(t_{n-1}, t_n)$, only the account value varies since all the benefit bases, death and life, remain constant. Thus, for $t \in (t_n^+, t_{n+1}^+]$, the annuity value $\Phi(t, A, G)$ solves the following linear PDE for each fixed value of the benefit base G

$$ \partial_t \Phi + \mathcal{L} \Phi = 0, \quad (2.15) $$

where the operator $\mathcal{L}$ is

$$ \mathcal{L} \Phi = \frac{1}{2} \sigma^2 A^2 \partial_{AA} \Phi + (r - \alpha^A) A \partial_A \Phi - \alpha^G G \partial_A \Phi - r \Phi. \quad (2.16) $$
Localization and boundary conditions

Equation (2.15) is originally posed on the domain $(t, A) \in [0, T] \times [0, \infty)$. For computational purposes, and because asset prices are finite and so is the account value, one needs to localize this domain to $[0, T] \times [0, A_{\max}]$ where $A_{\max}$ is large enough not to be attained by the account value during the lifetime of the annuity. Thus, we need to add complementary boundary conditions. We consider that we are between two anniversary dates $t_n^+$ and $t_{n+1}^-$ backwards.

— When $A = 0$, the policyholder has no longer the possibility to make any withdrawal from his account. However, if the IB election is possible, then the income period begins, given the policyholder is alive, and the death benefit is activated if he is dead at $t_{n+1}$. Since the account value is equal to zero, then the annuitization will be indexed on the benefit base. Therefore, we have:

$$ \Phi(t, 0, G) = e^{-r(t_{n+1}-t)} \left( p_{n+1} \frac{\ddot{a}_{t_{n+1}}^{gua}}{\ddot{a}_{t_{n+1}}^{act}} G + p_n q_{n+1} D(t_{n+1}^{-}, 0, G) \right). \quad (2.17) $$

— When $A = A_{\max}$, we consider retrieving all the cash more interesting than any other strategy if the policyholder is alive. If he dies, the death benefit will be activated. Therefore, the Dirichlet boundary condition for this case is:

$$ \Phi(t, A_{\max}, G) = e^{-r(t_{n+1}-t)} \left( p_{n+1}A_{\max} + p_n q_{n+1}A_{\max} \right) = e^{-r(t_{n+1}-t)} p_n A_{\max}. \quad (2.18) $$

Let us define the solution domains:

$$ \begin{aligned} \tilde{\Omega}_n &= [t_{n-1}^+, t_n^-] \times [0, A_{\max}] \\ \tilde{\Omega} &= \bigcup_{t_n} [t_{n-1}^+, t_n^-] \times [0, A_{\max}]. \end{aligned} $$

The pricing problem for the GMIB variable annuity combined with DB under the discrete withdrawal scenario is then achieved in $\tilde{\Omega}$ as follows: within each set $\tilde{\Omega}_n$, $n = 1, \dots, N-1$, the solution to the problem is the viscosity solution of a decoupled set of linear PDEs (2.12) with final condition (2.11) and boundary conditions (2.17)-(2.18) computed from the nonlinear algebraic Equation (2.14).

Construction of the scheme

Let $(A_0, A_1, \dots, A_J)$ be the equally spaced grid in the direction of the account value with $A_0 = 0$ and $A_J = A^{\max}$. Analogously $(G_0, \dots, G_K)$ is an equally spaced grid for the benefit base with $G_0 = 0$ and $G_K = G^{\max} = A^{\max}$. The spacial steps for both variables are considered to be equal. That is:

$$ \Delta A = \Delta G, $$

where $\Delta A = \frac{A_{\max} - A_0}{J}$ and $\Delta G = \frac{G_{\max} - G_0}{K}$.

Hence, $A_j = j\Delta A$ and $G_k = k\Delta G, \forall j, k$. The discrete time steps are denoted by $n\Delta t$ for $n = 1, \dots, N$ where $T = N\Delta t$. Since, in our analysis, we consider that events occur only at anniversary dates which are yearly, $\Delta t = 1$ and each time $t_n$ coincides with the discrete time step $t_n = n$.

The numerical procedure to solve the approximation in (2.15) is the standard finite difference approach. We employ the two-level implicit finite difference scheme to discretize the differential term $\mathcal{L}\Phi$ as given in (2.16). Let $\mathcal{L}_h\Phi_{j,k}^n$ denote the discrete value of the differential operator at the node $(A_j, G_k, t_n)$. The approximation is then given by:

$$ \mathcal{L}_h \Phi_{j,k}^n = \frac{\sigma^2}{2} A_j^2 \frac{\Phi_{j+1,k}^n - \Phi_{j,k}^n + \Phi_{j-1,k}^n}{\Delta A^2} + \{(r - \alpha^A)A_j - \alpha^G G_k\} \frac{\Phi_{j+1,k}^n - \Phi_{j-1,k}^n}{2\Delta A} - r \Phi_{j,k}^n $$
The general theta-scheme for solving Equation (2.15) is given by

$$ \frac{\Phi_{j,k}^{n+1} - \Phi_{j,k}^{n}}{\Delta t} = \theta \mathcal{L} \Phi_{j,k}^{n+1} + (1 - \theta) \mathcal{L}_h \Phi_{j,k}^{n} $$

where $\theta$ is a weighting factor, $0 < \theta \le 0$. For $\theta = 0$, this scheme is the explicit scheme, whereas $\theta = 1$ corresponds to the implicit one. The error of the previous cases is of $O(\Delta A^2, \Delta t)$. The explicit scheme has stability issues while the implicit scheme is absolutely stable. The most popular scheme for approximating the solution of the Black-Scholes equation is the Crank-Nicolson scheme obtained for $\theta = 1/2$. The latter is shown to be unconditionally stable and $O(\Delta A^2, \Delta t^2)$ convergent, see [68]. In particular, this scheme is used in [60] to solve the optimal pricing problem of the GMWB rider with rational behavior.

The discretization w.r.t the benefit base is not important here. However, since the PDE is solved backwards between $t_{n+1}^{-}$ and $t_{n}^{+}$, we need to divide this time period, i.e. the period between two consecutive withdrawal dates, into finer time steps for a good accuracy due to the finite difference approximation to the partial derivatives.

### Applying the jump condition

Recall that changes in the benefit base only occur at withdrawal dates. After withdrawing the amount $\gamma_n$ at time $t_n$, the account value changes from $A_{t_n^{-}}$ to $A_{t_n^{+}} = h^A(A, G, \gamma_n)$, and the benefit base drops from $G_{t_n^{-}}$ to $G_{t_n^{+}} = h^G(A, G, \gamma_n)$. The jump condition of $\Phi(t_n, A, G)$ across $t_n$ is given by

$$ \Phi(t_n^{-}, A, G) = \max_{0 \le \gamma_n \le A} \left( \Phi(t_n^{+}, h^A(A, G, \gamma_n), h^G(A, G, \gamma_n)) + p_n \bar{f}(A, G, \gamma_n) + p_{n-1} q_n D(t_n^{-}, A, G) \right) \quad (2.19) $$

For the optimal strategy, the withdrawal amount $\gamma_n$ is chosen under the restriction $0 \le \gamma_n \le A$ to maximize the value of $\Phi(t_n^{-}, A, G)$ in Equation (2.19).

The application of the jump condition decreases the account value and benefit base. For each $G_j$, a continuous solution from PDE (2.15) is associated. We can restrict the possible values for the withdrawal amount to multiples of $\Delta A$. This implies, for a given account value $A_j$ at time $t_n^{-}$, the withdrawal amount $\gamma$ takes $j$ possible values: $\gamma = A_j - A_i$, $l = 1, 2, \dots, j$. However, numerical tests showed that a finer grid is preferable for the withdrawal amount. Therefore, it is not guaranteed that the account value, nor the benefit base after the withdrawal, $A_{t_n^{+}}$ and $G_{t_n^{+}}$, fall within their respective grid nodes. To solve this issue, a two-dimensional interpolation is required. In this work we adopted a bi-linear interpolation.

Suppose the jump condition requires the value $\Phi(., A, G)$ at the point $(A, G)$ located inside a grid $A_i \le A \le A_{i+1}$ and $G_j \le G \le G_{j+1}$, then the interpolation is performed as the following:

$$ \Phi(., A, G) \approx \frac{G_{i+1} - G}{G_{i+1} - G_i} \Phi(., A, G_j) + \frac{G - G_i}{G_{i+1} - G_i} \Phi(., A, G_j), \quad (2.20) $$

where:

$$ \Phi(., A, G_j) \approx \frac{A_{i+1} - A}{A_{i+1} - A_i} \Phi(., A_i, G_j) + \frac{A - A_i}{A_{i+1} - A_i} \Phi(., A_{i+1}, G_j), $$

$$ \Phi(., A, G_j) \approx \frac{A_{i+1} - A}{A_{i+1} - A_i} \Phi(., A_i, G_{j+1}) + \frac{A - A_i}{A_{i+1} - A_i} \Phi(., A_{i+1}, G_{j+1}). $$

At last, the jump condition is achieved through combining (2.19) and (2.20) to find the optimal withdrawal and maximize the function $\Phi$.
## Similarity and dimension reduction

An important feature of the contract value is that it exhibits good scaling properties in the Black-Scholes case. We can easily verify that the solution $\Phi(t, A, G)$ of PDE (2.15) with boundary conditions (2.17) and event conditions (2.18) verifies

$$ \Phi(t, \xi A, \xi G) = \xi \Phi(t, A, G) $$

for any scalar $\xi > 0$. Therefore, choosing $\xi = 1/G$ we obtain

$$ \Phi(t, A, G) = G\Phi(t, \frac{A}{G}, 1) = G\Phi(t, \tilde{A}), $$

where $\tilde{A} = \frac{A}{G}$. It means that we need only solve the corresponding equations for the one-dimensional function $\phi$ defined in the following:

— Between two consecutive withdrawal dates $(t_{n-1}, t_n)$ $\phi$ follows the PDE

$$ \partial_t \phi + \frac{1}{2} \sigma^2 \tilde{A}^2 \partial_{\tilde{A}\tilde{A}} \phi + (r - \alpha^A) \tilde{A} \partial_{\tilde{A}} \phi - \alpha^G \partial_{\tilde{A}} \phi - r \phi = 0 \quad (2.21) $$

— At the anniversary dates $t_n$, the jump condition is explicitly expressed as the following:

$$ \phi(t_n^-, \tilde{A}_{t_n^-, \tilde{Y}_n}) = \max_{\tilde{Y}} \left( h_1(\tilde{Y})\phi(t_n^+, h_2(\tilde{A}_{t_n^+, \tilde{Y}_n})) + p_n \tilde{Y}_n + p_{n-1} q_n D(t_n^+, \tilde{A}_{t_n^+, 1}) \right) $$

where the functions $h_1$ and $h_2$ are a reduced version of the account value and benefit base evolution at the anniversary dates, and are defined according to different cases:

1. Roll-up only case:

$$ h_1(\tilde{Y}) = \begin{cases} 1 + \eta - \tilde{Y}_n & \text{if } \tilde{Y}_n \le \eta \\ 1 - \frac{\tilde{Y}_n - \eta}{\tilde{A}} & \text{if } \tilde{A} \ge \tilde{Y}_n > \eta \end{cases} $$

$$ h_2(\tilde{Y}) = \begin{cases} \frac{\tilde{A} - \tilde{Y}_n}{1 + \eta - \tilde{Y}_n} & \text{if } \tilde{Y}_n \le \eta \\ \frac{\tilde{A} - \tilde{Y}_n}{1 - \frac{\tilde{Y}_n - \eta}{\tilde{A}}} & \text{if } \tilde{A} \ge \tilde{Y}_n > \eta \end{cases} $$

2. Ratchet only case:

$$ h_1(\tilde{Y}) = \max(\tilde{A} - \tilde{Y}_n, 1 - \frac{\tilde{Y}_n}{\tilde{A}}) $$

$$ h_2(\tilde{Y}) = \min(1, \tilde{A}) $$

3. Reset case (Roll-up + ratchet):

$$ h_1(\tilde{Y}) = \begin{cases} \max(\tilde{A} - \tilde{Y}_n, 1 + \eta - \tilde{Y}_n) & \text{if } \tilde{Y}_n \le \eta \\ \max(\tilde{A} - \tilde{Y}_n, 1 - \frac{\tilde{Y}_n - \eta}{\tilde{A}}) & \text{if } \tilde{A} \ge \tilde{Y}_n > \eta \end{cases} $$

$$ h_2(\tilde{Y}) = \begin{cases} \min(1, \frac{\tilde{A} - \tilde{Y}_n}{1 + \eta - \tilde{Y}_n}) & \text{if } \tilde{Y}_n \le \eta \\ \min(1, \frac{\tilde{A} - \tilde{Y}_n}{1 - \frac{\tilde{Y}_n - \eta}{\tilde{A}}}) & \text{if } \tilde{A} \ge \tilde{Y}_n > \eta \end{cases} $$

It can be easily verified that PDE (2.21) does not depend on the benefit base since the latter is constant between two consecutive withdrawal dates $t_{n-1}^+$ and $t_n^-$. Therefore, the resolution of the two-dimensional problem can be reduced into a one-dimensional problem. This is particularly useful when adding more stochastic variables like stochastic volatility and/or interest rates.
## 2.5 Results

The main goal of this study is to assess the behavior risk of a given GMIB product, in case policyholders maximize their expected cash flows. Through optimal withdrawal amounts, or IB election, the insurer is concerned that his product may lead to "undesirable" policyholders behaviors. Given a state of the universe in a future time, and a set of up-front fixed hypotheses (management and guarantee fees, interest rates and volatility), the optimal framework allows us to predict these behaviors through the stochastic control problem detailed in previous sections. As of the product, we consider Retirement Cornerstone©commercialized by AXA U.S. The product hypotheses usually change to account for a new financial environment or customers needs.

We choose two close variations of the product launched in the last decade, which we call them respectively, Product A and Product B. These are two GMIB products to which a death benefit (DB) can be added, i.e. Product A-DB and Product B-DB. The parameters values are given in Table 2.5 :

<table><caption>TABLEAU 2.5 – Model parameters</caption><thead><tr><td>Parameters</td><td>Product A</td><td>Product B</td></tr></thead><tbody><tr><td>Policyholders initial age x<sub>0</sub></td><td>60</td><td>60</td></tr><tr><td>First date for IB election</td><td>10th</td><td>10th</td></tr><tr><td>Last age for IB election</td><td>85th</td><td>95th</td></tr><tr><td>Last age for DB if any</td><td>95</td><td>95</td></tr><tr><td>Interest rate r</td><td>2% and 4%</td><td>2% and 4%</td></tr><tr><td>Volatility σ</td><td>20%</td><td>20%</td></tr><tr><td>Roll-up rate η</td><td>6%</td><td>r + 1%</td></tr><tr><td>Initial premium A<sub>0</sub></td><td>$100,000</td><td>$100,000</td></tr><tr><td>Total fees ᾱ = α<sup>A</sup> + α<sup>G</sup></td><td>3.5%</td><td>3.5%</td></tr></tbody></table>

We compare these variations for roll-up and reset, with and without death benefit. We conduct the experiments to illustrate the following:

- Policyholders rational behavior for Products A, B, A-DB and B-DB based on the two dimension approach, giving the withdrawal surface as a function of the account value A and benefit base G, for two different interest rate values, in four different periods in time. We also give withdrawals as a function of the moneyness A/G based on the dimension reduction approach.

- Policyholders rational behavior for different values of the volatility σ = {10%, 20%, 30%, 40%}.

- The contract initial value as a function of interest rates for Product A-DB.

- The contract initial value as a function of the total fees.

### Overview of the policyholder behavior

#### Roll-up only case

First, in Figure 2.2, we present the withdrawal amount surface as a function of the account value A and benefit base G, for Products A and A-DB with roll-up only for fixed times $t = 3, 13, 23$ and $r = 2\%, 4\%$.
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>FIGURE 2.2 – Policyholder optimal withdrawal amount as a function of the account value A and benefit base G for Product A and A-DB</figcaption></figure>
As a first remark, the guaranteed account must be higher than the account value for the policyholder to stay in the contract. Depending on the moneyness, he or she can choose not to withdraw, or withdraw the guaranteed amount. The IB election (using the guaranteed account) is left at maturity or for very small account values. These finding are actually confirmed in Figure 2.3 where we give the withdrawal strategy as a function of time $t$ and moneyness A/G based on the dimension reduction scheme.

Product A

<figure><img src="image_24.png" /><figcaption>FIGURE 2.3 – The withdrawal strategy $\tilde{Y} = \frac{Y}{G}$ as a function of time $t$ and moneyness A/G</figcaption></figure>

Based on Figures 2.2-2.3, we can point out the following remarks :

- There is a wider range of guaranteed withdrawals for $r = 4\%$ compared to $r = 2\%$. It tells us that for a roll-up rate as high as 6%, the policyholder tends to wait for the benefit base to increase at this rate in a low interest rates environment instead of draining the subaccount and guaranteed account.

- The death benefit increases the expected cash flows in the future, which is also a motif for the policyholder to wait.

- The IB election indexed on the guaranteed account is only expected to happen in the absence of a death benefit. Even in such case, it only takes place closer to the maturity of the product and for small account values. Actually, the ratio $\frac{\dot{g}^{act}}{\dot{g}ua}$ is not very favorable for the insured, and he or she would rather start withdrawing the guaranteed amount few years earlier.
— Actually, fees are quite high so the account value usually drops quickly. This restricts the analysis to A/G relatively small.

Note that reducing the dimension allows to increase the speed of calculations.

Product A was launched in a period where interest rates were around 4% which justifies the choice of a 6% roll-up rate. Given the behavior expressed in the Figures 2.3, the actual interest rate level may seem quite high and the product more interesting for the insured than the insurer.

Later, Product B was launched with reconsidered assumptions. The roll-up rate becomes indexed on the short term interest rates with a spread of 1%. The insurance company attracted the costumer by setting a longer limiting age to annuitize (until the policyholder's 95th anniversary). In Figure 2.4 we give the policyholder behavior in time as a function of the moneyness.

Product B

<figure><img src="image_25.png" /><figcaption>FIGURE 2.4 – The withdrawal strategy as a function of time *t* and moneyness A/G</figcaption></figure>

While the behavior remains very close, we can however notice that in the absence of the death benefit, partial withdrawals disappear, the likelihood of the IB election is higher for low interest rate, and left until the last years for higher interest rate. Moreover, the likelihood of recovering the account value
is higher than in Product A. In the presence of a death benefit withdrawing the guarantee becomes interesting around the 20th anniversary of the contract, while an IB election is exercised at maturity.

## Reset case

Ratchets allow for the benefit base to be set at the account value level when the latter is higher than the previous benefit base level. Combined with the roll-up, we have the reset which is a very attractive feature for policyholders who are interested in the stochastic performance of stock markets, but at the same time want to have a guaranteed minimum performance. We present in Figures 2.5 and 2.6 the results related to products A and B with and without death benefit for the previous interest rates and roll-up values. In this case, the insured sticks with guaranteed or zero withdrawals for most time, and

### Product A

<figure><img src="image_26.png" /><figcaption>FIGURE 2.5 – The withdrawal strategy as a function of time $t$ and moneyness $A/G$</figcaption></figure>

tends to elect the income benefit in the last anniversary date if the account value is low, and recover it otherwise. The reset is very costly for the insurance company, however, fees are quite high for this product and the ratchet takes place only at early dates since the account value is brought down by the fees rate.

In what follows, we will focus on the roll-up only case in an attempt to analyze the impact of some key parameters in the pricing and expected policyholders rational behavior for these products.
Product B

<figure><img src="image_27.png" /><figcaption>FIGURE 2.6 – The withdrawal strategy as a function of time *t* and moneyness A/G</figcaption></figure>

The impact of volatility

The volatility level assumption is very important for variable annuities in general and the GMIB pro-
duct in particular. Based on Product A for the hypothesis used above, we compare two levels of vo-
latility (which we can compare to the 20% volatility case given in Figure 2.3). We can see in Figures
2.7 and 2.8 that the lower the volatility, the earlier guaranteed withdrawals start. Moreover, the lapsing
likelihood is also higher. This means that, the more uncertain are markets, the more the policyholder
tends to withdraw money from his account. On the other hand, the IB election does not seem to be
affected.

Roll-up rate and fees

There is a trade-off between roll-up rate and fees. The roll-up is the mechanism that ensures the policyholder a minimum return, which can be higher than the money market. However, to be able to provide interesting roll-up rates, insurance companies need to be hedged from uncertain interest rates. In this sense, they use for example swaps. Therefore, they need to collect fees that at least allow for a fair pricing for the contract, i.e. such that the paid cash flows equal the premium. On the other hand, high fees can have a perverse effect. By decreasing the subaccount value, especially in periods of low performance, present collected fees reduce future potential ones. On the long run, combined with
$\sigma = 10\%$

<figure><img src="image_28.png" /><figcaption>FIGURE 2.7 – The withdrawal strategy as a function of time $t$ and moneyness $A/G$ for $\sigma = 10\%$</figcaption></figure>

guaranteed withdrawals, the income benefit can be elected by bringing the account value to zero. Moreover, when the account value falls to zero, the insurance company can no longer collect fees and starts to pay the guarantee.

In Figures 2.9, we compare the value of the contract at inception for different parameters of Product B as a function of total fees. The fair price corresponds to $\phi(0, \bar{A}_0) = 1$. We see that the contract is underpriced with death benefit. The fair fees would be as high as 7%. Without the death benefit, they are around 3% which is close to the rates applied by the insurance company. In Figure 2.10, we conduct a similar test by varying the roll-up rate for Products A, A-DB, B and B-DB for $r = 2\%$. We notice that $\phi(0, \bar{A}_0) = 1$ corresponds to a roll-up rate that is close to the interest rate except one of the products. Indeed, Products A, A-DB and B-DB are under-priced for the features they provide. On the other hand, the insurance company was conservative in the roll-up rate assumption for Product B which allows it to be profitable even for the worst case scenario. Of course, insurance companies do not expect (and hope not) that all policyholders follow an optimal behavior. However, prudent hypotheses can prevent from important losses. Including a proportion of policyholders that are likely to behave optimally is one of the solutions. Note that the GMIB product is less risky than the GMWB in that the annuity factor is defined with conservative assumptions.
$\sigma = 30\%$

<figure><img src="image_29.png" /><figcaption>FIGURE 2.8 – The withdrawal strategy as a function of time $t$ and moneyness $A/G$ for $\sigma = 30\%$</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>FIGURE 2.9 – The contract value at inception as a function of the total fees $\bar{\alpha}$ for Product A with and without DB for $r = 2\%$, $4\%$</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>FIGURE 2.10 – The contract value at inception as a function of the roll-up rate $\eta$ for Product A and B with and without DB for $r = 2\%$</figcaption></figure>
## 2.6 Conclusion

In this work, we analyzed the optimal behavior of a policyholder entering a GMIB contract combined with a death benefit guarantee. The solution is based on a optimal stochastic control framework in a Black-Scholes framework, and solved numerically using recursive dynamic programming techniques. Considering only the account value varies between two anniversary dates, we used finite differences methods and a linear search for the optimal withdrawal to maximize the expectation of discounted future cash flows. Such calculations give an optimal withdrawal function that depends on time, account value and benefit base. Taking advantage of the good scaling properties provided by the contract payoff and the asset price, we are able to reduce the dimensionality of the problem to time and moneyness, making calculations faster and results interpretation easier.

The policyholder's optimal behavior corresponds to the worst case scenario for insurance companies. Therefore, even though insurers are not expected to behave optimally, a good insight of how they may act in case they do, given a market environment, can allow insurers to be more effective in pricing and hedging their products. We find that the optimality consisted mainly in four choices; zero withdrawals, guaranteed withdrawals, lapse and IB election. We presented these results for two different products before analyzing the impact of some of the contract key parameters.

Finally, these results can be used as a guide to practitioners in the design of new products where a particular client behavior is desired. The model can be used to compute the fair withdrawal fee with hedging purposes. In particular, we find that these products are under-priced in case of an optimal behaviors, as it was already mentioned by [130] in the case of GMWBs. We believe that understanding policyholder behavior is a critical concern to the future of life insurance business, and due to its importance, further research is required.
# Chapitre 3

## iCPPI and Gap Risk

***Abstract***— Individualized Constant Proportion Portfolio Insurance (iCPPI) products are attractive alternatives to traditional unit linked products offering a guaranteed minimum return, such as variable annuities. They offer high potential returns whilst limiting the downside risk by implementing a dynamic allocation strategy between risky and risk-free assets tailored to the risk attitude of the beneficiary. But performance evaluation of iCPPI products should not rely on the unrealistic assumptions of continuous market price variations and continuous rebalancing of asset allocations. We adopt a more general and realistic price jump model and examine several dynamic strategies as well as gap put options to mitigate the risk that the value of the product falls below the guaranteed minimum.

**Keywords**: CPPI; dynamic multiplier; jump processes and gap risk; vanilla and gap options.

### 3.1 Introduction

Increased market volatility and falling interest rates triggered by the 2008-2009 financial crisis reduced the performance of traditional long-term investment products, increased their risks and, where applicable, their capital requirements. In this context the new iCPPI products provide an attractive alternative to many traditional long-term investment products offering a guaranteed minimum return, such as variable annuities, for several reasons : lower exposure to uncertain volatilities and extreme market price movements, lower costs, and lower regulatory capital requirements, to name a few.

Already, with rising life expectancy, current provisions for retirement may not be sufficient for many people to secure acceptable life standards after retirement. To achieve sufficiently high investment returns, together with low risks over the long term, funds should remain invested in stocks and other risky assets as well as in the safer bonds over an extended period well into retirement. The design of long-term investment products should also reflect the requirements and risk attitudes of individual investors.

Constant Proportion Portfolio Insurance (CPPI) is the name given to an investment strategy that provides a minimum guaranteed return, the “floor” (usually defined as the discounted value of a final capital guarantee) and aims to maintain, at all times, an exposure to a risky asset equal to a constant multiple of the “cushion” defined as the excess value of the fund above the floor. The final capital guarantee and the multiplier are chosen to satisfy the risk attitude of the investor.

Assuming the risky asset follows a geometric Brownian process, constant rate for the risk-free asset, and continuous relancing with no transaction costs, the CPPI payoff is optimal for an investor with
a risk tolerance coefficient varying linearly with wealth, see [29], [125], [140]. Specifically, the CPPI payoff is equal to the floor plus a cushion value which is proportional to asset price power the multiplier. The floor and multiplier are chosen according to the two parameters of a HARA utility function to maximize the expected utility of the investor. Additional advantages offered by CPPI strategies over more traditional investments with minimum guaranteed returns are: price transparency, open time-horizon, no early redemption penalty, wide range of alternative investments for the risky asset, and flexibility to add other guarantees such as ratchets (see II.1.4).

iCPPI is a CPPI strategy adapted to evolving individual needs and market conditions. The floor and multiplier are modified accordingly. Thus iCPPI may combine most of the advantages of CPPI with the need for flexibility and enhanced risk management.

ICPPI providers (typically, an insurance company) face many challenges in the implementation of the dynamic strategy that replicates the guaranteed payoff. Adjusting the risky asset/ risk-free asset allocation can take place at discrete times only, there are transaction costs, and risky asset prices may jump. Thus, there is a difference between the theoretical value of an iCPPI strategy under hypothetical assumptions, and reality. In particular, there is a positive probability for the value of the fund to fall below the guaranteed floor. We call such shortfall the gap risk.

The analysis of the gap risk has often been limited to simple conditions to preserve analytical tractability:

- Unrealistic modeling of the risky asset price market including continuous price dynamics, zero-cost trading and unlimited liquidity.

- Simple parameterization of the CPPI strategy such as constant capital guarantee and multiplier.

- Simplistic rebalancing strategies such as constant frequency.

As a result, the iCPPI offers a mechanism that takes advantage of the specific advantages of both stocks and bonds, while complying with growing needs of flexibility as experienced by policyholders.

However, the implementation of iCPPIs at insurance companies levels suffers from a number of operational constraints on the asset management: the rebalancing occurs through regular albeit discontinuous (at most daily) checks between the insurance company and a bank; depending on the design of the iCPPI and the discontinuous rebalancing frequency, the magnitude of the earnings at extreme risk may require the externalization of the gap risk management to the bank. As a result the main issue experienced by the insurance company remains to minimize the downside risk and keep control of the gap risk, which involves three main challenges: This article extends previous analyses of the gap risk by introducing:

- Price jump dynamics

- A dynamically adjusted multiplier

- Advanced rebalancing strategies, vanilla and gap put options to mitigate the gap risk

### 3.1.1 A brief review of the literature

CPPI as a mechanism falls within portfolio management techniques which ensure a lower bound on the portfolio value at a given maturity. In theory, one can have protection against unfavorable market scenarios for some asset by investing in a put option with strike equal to the desired lower bound. This is known as the option based portfolio insurance (OBPI) and was introduced by [116] and [41]. However, the put option needed to perform this protection may not be available in the market, for example
if the investment horizon is long. And even though one can replicate the payoff of the put option by trading the asset and the cash, such replication is costly and imperfect.

The CPPI was introduced by [140] and [29]. Its wide use in the financial industry, see [138], brought a lot of attention. Typical buyers are large individual investors and institutional investors such as pension funds. The main topics the literature on the CPPI covers are: limiting its risk, developing hedging techniques to cover the remaining risk, and the behavior of the CPPI.

[26] defined an upper bound for the multiplier $m$ such that the investment in the risky portfolio is maximized under the gap risk must stay under a certain limit. [45] extended the calculation of the VaR and GVar of the CPPI portfolio. [144] and [57] on the other hand, contributed in extending the CPPI in order to build a protection against the small but existing gap risk. The price and size of such protection is model-dependent and will depend on the probability of hitting the floor. Finally, [26] and [84] study in detail the behavior of the CPPI strategy under specific conditions for the underlying asset. [26] for example considered the case where the underlying risky asset of the CPPI fund is an index or a basket of indexes. They used a Multivariate Variance Gamma (MVG) model for a series of correlated spreads to price the CPPI.

### 3.1.2 Review of CPPI mechanism basics

Consider at time $t$ a risky asset (e.g., a share) with price $S_t$ and a risk-free asset (e.g., a Treasury bond) with price $B_t$ returning a constant rate $r$. The CPPI fund is invested into these two assets so that part of its value, called the "floor" $F_t$, is guaranteed whilst the excess value above the floor, called the "cushion" $C_t = V_t - F_t$, remains exposed to the risky asset price fluctuations. At any time, the exposure to the risky asset, $e_t$, is kept at a constant multiple, $m$, of the cushion, that is:

$$e_t = mC_t$$

The rest of the value of the fund is invested (or, if negative, borrowed) at the risk-free rate (Note that the exposure $e_t$ may be acquired at no cost if using an off-balance sheet instrument such as a future, which may be advantageous because of liquidity and low transaction costs). The floor is often chosen to increase over time at the risk-free rate (it could not be made to increase faster indefinitely), that is:

$$F_t = F_0e^{rt}$$

In theory, when the risky asset price follows a geometric Brownian motion, and with continuous, zero-cost rebalancing (Black-Scholes conditions), the value of the cushion is path independent and proportional to $S_t^m$. In other words, it is the value of a power option. It is convex when $m > 1$ (like a long call option), linear when $m = 1$, and concave when $m < 1$, like a short put option. But unlike standard call and put options there is no need to fix an expiry date, a CPPI strategy is open-ended. Under the above assumptions, the value of the cushion would never fall to zero; in practice, if it does fall to zero or below zero (e.g., because of a price jump or of discrete rebalancing), the entire fund is monetized, i.e., is entirely invested in the risk-free asset, and the product provider must make up the shortfall to deliver the floor value. In practice there may also be other constraints such as no borrowing or additional features such as ratcheting up the floor. In those cases, the path independency and open-endedness of the product are lost and the payoff profile becomes more complex.
## 3.2 Methodolody

### 3.2.1 CPPI in theory and practice

#### Continuous-time framework

The risky asset S is defined by the diffusion equation $dS_t = \mu S_t dt + \sigma S_t dW_t$ where W is a standard Brownian motion. The previous hypothesis for the risk-free asset are kept. In such context, and assuming continuous time CPPI, the cushion C is log-normally distributed with drift $m(\mu - r) + r$ and volatility $m\sigma$:

$$C_t = C_0 \exp\left(\left(m(\mu - r) + r - \frac{m^2\sigma^2}{2}\right)t + m\sigma W_t\right)$$

and the portfolio value V has the path independent expression:

$$V_t = F_t + (V_0 - F_0) \exp\left(\left(m(\mu - r) + r - \frac{m^2\sigma^2}{2}\right)t + m\sigma W_t\right)$$

However, such assumptions are unrealistic and not consistent with market practice. To remedy these unrealistic hypothesis, two alternatives are studied: modeling in a discrete-time framework and in a Lévy framework.

#### Discrete-time CPPI

In practice the CPPI is rebalanced in discrete time, where the shortfall probability is no longer equal to 0, which implies to monetize more often.

A sequence of equidistant refinements of the interval $[0, T]$ is defined:

$$\Theta = \{t_0 = 0 < \dots < t_{N-1} < t_N = T\}$$

where $t_{k+1}^N - t_k^N = \frac{T}{N}$ for $k = 0, \dots, N-1$. The number of shares is constant on the intervals $]t_i, t_{i+1}]$. Let $t_s := \min\{t_k \in \Theta | V_{t_k} - F_{t_k} \le 0\}$. The first time the portfolio value touches the floor. The discrete-time cushion follows the equation:

$$C_{t_{k+1}} = e^{r(t_{k+1}-\min\{t_s, t_{k+1}\})} (V_{t_0}^{\Theta} - F_{t_0}) \prod_{i=1}^{\min\{s, k+1\}} \left( m \frac{S_{t_i}}{S_{t_{i-1}}} - (m-1)e^{r\frac{T}{N}} \right),$$

or recursively:

$$C_{t_{k+1}} = \begin{cases} C_{t_k} \left( m \frac{S_{t_{k+1}}}{S_{t_k}} - (m-1)e^{r\frac{T}{N}} \right) & \text{if } C_{t_k} > 0, \\ C_{t_k} e^{r\frac{T}{N}} & \text{if } C_{t_k} \le 0. \end{cases} \quad (3.1)$$

$V_{t_k}$ is given through the relation $V_{t_k} = C_{t_k} + F_{t_k}$.

To comply with the CPPI algorithm and respect practical constraints, the number of shares of the risky and safe assets ($\alpha$ and $\beta$) are as follows:

$$- \alpha_{t_k} = \min\left(\max\left(\frac{mC_{t_k}}{S_{t_k}}, 0\right), \frac{V_{t_k}}{S_{t_k}}\right).$$

$$- \beta_{t_k} = \frac{V_{t_k} - \alpha_{t_k} S_{t_k}}{B_{t_k}}.$$

When adding transaction costs, they are taken as a proportion of the change in the risky exposure, i.e $\propto (\alpha_{t_k} - \alpha_{t_{k-1}}) \times S_{t_k}$. So at time $t_k$, the number of shares of the risky asset will be reduced to:

$$\tilde{\alpha}_{t_k} = \alpha_{t_k} - |\alpha_{t_k} - \alpha_{t_{k-1}}| \times \text{nb of bps}$$
The CPPI capital guarantee is ensured, as long as the bond floor is not breached through enabling to fully invest the portfolio into the non risky assets. The probability of breaching the floor is defined as the probability that the portfolio value falls below the floor, i.e $P^{BF} := \mathbb{P}(V_T \le G) = \mathbb{P}(\exists t \in [0, T] : V_t \le F_t)$. The local shortfall probability is the conditional probability defined as : $P_{t_i, t_{i+1}}^{LBF} = \mathbb{P}(V_{t_{i+1}} \le F_{t_{i+1}} | V_{t_i} > F_{t_i})$. The two are related as follows $P^{BF} = 1 - \prod_{i=1}^{i=N} (1 - P_{t_i, t_{i+1}}^{LBF})$.

This probability which was equal to zero in the continuous Black-Scholes model, is now greater than zero. Assuming the portfolio did not breach the floor up to $t_k$, the probability of breaching the floor at $t_{k+1}$, is that of a downside jump in the risky asset of more than about $1/m$. Its mathematical expression is :

$$P_{t_i, t_{i+1}}^{LSF} = \mathbb{P}\left(\frac{S_{t_{i+1}}}{S_{t_i}} \le \frac{m-1}{m} e^{\frac{T}{N}}\right),$$

where the evolution of the risk-free part with rate $r$ is taken into account.

The backtesting is based on the period Q1-2006 to Q4-2010 on S&P500 index. Simulating paths (N=10,000) in the Black & Scholes model is made using the 3-month realized volatility based on the standard deviation (see Figure 3.1), a constant asset return $\mu = 8\%$. The rate of the risk-free asset is $r = 4\%$. Three rebalancing frequencies are being compared regarding the distribution of the final portfolio value (daily, weekly and monthly), with the following assumptions :

- Initial investment/Guarantee : $100, and $100

- Duration : 5 years

- Transaction costs : 10 bps

<table><caption>TABLEAU 3.1 – Final value metrics : Buy & Hold strategy vs CPPI with m = 3 vs CPPI with m = 6</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Buy & Hold Strategy</th><th colspan="3">CPPI with m = 3</th><th colspan="3">CPPI with m = 6</th></tr><tr><th>Daily</th><th>Weekly</th><th>Monthly</th><th>Daily</th><th>Weekly</th><th>Monthly</th></tr></thead><tbody><tr><th>Mean</th><td>126.97</td><td>123.31</td><td>122.39</td><td>119.75</td><td>124.10</td><td>124.87</td><td>125.01</td></tr><tr><th>Std-Dev</th><td>7.18</td><td>31.58</td><td>32.66</td><td>36.86</td><td>42.62</td><td>43.88</td><td>48.10</td></tr><tr><th>95% quantile</th><td>116.90</td><td>100.48</td><td>99.98</td><td>97.01</td><td>99.99</td><td>99.13</td><td>89.69</td></tr><tr><th>99.5% quantile</th><td>113.42</td><td>100.02</td><td>99.88</td><td>91.47</td><td>99.98</td><td>95.20</td><td>74.28</td></tr><tr><th>5% quantile</th><td>140.21</td><td>194.37</td><td>195.23</td><td>197.94</td><td>216.51</td><td>218.50</td><td>225.46</td></tr><tr><th>0.5% quantile</th><td>150.63</td><td>266.47</td><td>284.07</td><td>282.58</td><td>291.49</td><td>293.75</td><td>311.46</td></tr><tr><th>Rebalancing cost</th><td>0.01</td><td>0.91</td><td>0.44</td><td>0.26</td><td>0.78</td><td>0.46</td><td>0.31</td></tr><tr><th>P<sup>BF</sup></th><td>0</td><td>0.0018</td><td>0.0947</td><td>0.5289</td><td>0.2016</td><td>0.5730</td><td>0.6555</td></tr></tbody></table>

The CPPI strategy under daily rebalancing performs better against a bear market than the weekly and monthly ones due to its reactiveness to decrease the risky exposure whenever needed. With such frequency, the guarantee is almost ensured; the less frequent we rebalance the more we are exposed to breaching the floor (as illustrated by fatter left tails (see Figure 3.2 bottom, right). The backtesting (Figure 3.2 top) and Table 3.5 illustrate the following remarks :

- In periods of mild market conditions, transaction costs negatively affect the performance of a daily rebalancing, although not to a significant extent.

- During a market crash, the three strategies monetize, with the daily rebalancing having less losses than the two others.

- The empirical probability of breaching the floor decreases when the rebalancing frequency increases.
— The cost of rebalancing increases with the frequency and with the multiplier. However, in our results, the cost of daily rebalancing for $m = 6$ is lower than the one with $m = 3$. This is explained by the fact that such a high multiplier allows for a total risky exposure and thus no rebalancing reducing the cost.

When comparing different strategies (Buy & Hold, CPPI with $m = 3$ and $m = 6$), we have the following results :

— The Buy & Hold strategy has higher expectation and lower standard deviation (table 3.5). This is mainly due to the low exposure to the risky asset. Its performance is highly correlated to the non-risky return (chosen to be 4%).

— The 5% and 0.5% quantiles show that the CPPI with $m = 6$ has a larger right tail and thus, performs better than the two others in bullish market. This remark is also illustrated in Figure 3.8.

Daily rebalancing almost prevents the bond floor from being breached, which ensures the capital guarantee at maturity. However, constant volatility and log-normal distribution modeling are not consistent with empirically observed jumps during extreme market moves likely to breach the bond floor. In order to relax these unrealistic assumptions, jumps are thus added through Lévy processes as developed in the next section.

## Adding jumps

We assume that the process of the risky asset follows a Lévy process :

$$ \frac{dS_t}{S_t} = dZ_t, $$

where $Z$ is a Lévy process. The risk-free asset $F_t$ is still deterministic.

Let $\tau = \inf\{t: V_t \le B_t\}$ the time where the portfolio value is fully invested in the risk-free asset. Until $\tau$ the actualized cushion ($C_t^* = \frac{C_t}{F_t}$) is as follows : $C_t^* = C_0^* \mathcal{E}(mL)_t$, where $\mathcal{E}$ denoting the stochastic exponential :

$$ \mathcal{E}(Z)_t = Z_0 e^{Z_t - \frac{1}{2} \langle Z \rangle_t} \prod_{s \le t, \Delta Z_s \ne 0} (1 + \Delta Z_s) e^{-\Delta Z_s}, $$

which gives us the portfolio value :

$$ V_t = \begin{cases} V_t \left\{ 1 + \left( \frac{V_0}{F_0} - 1 \right) \mathcal{E}(mL)_t \right\} & \forall t \le \tau, \\ V_\tau e^{r(t-\tau)} & \text{if } t > \tau. \end{cases} $$

The probability of breaching the floor can be expressed as :

$$ P^{\text{BF}} = \mathbb{P}(\exists t \in [0, T], V_t \le B_t) = 1 - \mathbb{P}\left(\forall t, \Delta L_t < \frac{1}{m}\right), $$

$$ P^{\text{BF}} = 1 - \exp\left(-T \int_{-\infty}^{-1/m} v(dx)\right), $$

which is illustrated by the fact that the number of downside jumps of size more than $\frac{1}{m}$ follows a Poisson distribution with intensity $T\nu(-\infty, -1/m)$.

For computation tractability, we choose the double exponential Kou model, see [112]). Under the risk neutral probability, the risky asset is modeled as follows :

$$ \frac{dS_t}{S_t} = \mu dt + \sigma dW + d\left(\sum_{i=1}^{N_t} e^{Y_i} - 1\right), $$
where W is a standard brownian motion, N is a poisson process with rate $\lambda$, the constants $\mu$ and $\sigma > 0$ are drift and volatility of the diffusion part and the jump sizes $\{Y_1, Y_2, ...\}$ are i.d.d random variables with a common asymmetric double exponential distribution of density:

$$f_Y(y) = (1 - p)\eta^+ e^{-\eta^+ y} 1_{y \ge 0} + p\eta^- e^{-\eta^- y} 1_{y < 0},$$

where $\eta^+$ is intensity of positive jumps while $\eta^-$ and $p$ are the intensity of negative jumps and the probability of their occurrence.

Under this jump model, and assuming a continuous rebalancing frequency, the probability of breaching the floor takes the following form:

$$P^{\text{BF}} = 1 - \exp\left(-Tp\lambda\left(1 - \frac{1}{m}\right)^{\frac{1}{\eta}}\right).$$

We substitute the Black-Scholes framework with the Kou model which we calibrate on implied volatility smile (between 2006 and 2011 for a 1-month implied volatility on a weekly basis ¹). We carry out the calibration by minimizing the quadratic error:

$$\sum_{i=1}^{9} \left( C_{t_i}(T, K_i)^{\text{Market}} - C_{t_i}^{\text{Kou}}(T, K_i, \sigma, p, \eta^+, \eta^-, \lambda) \right)^2,$$

where T is 1-month maturity, nine strikes $K_i$ from 80 to 110 and $(p, \eta^+, \eta^-, \lambda, \sigma)$ are the jump parameters (more details about. We give different statistics for these parameters in Table 3.2).

In order to avoid instability in parameters, we chose several starting points and set boundary condi-

<table><caption>TABLEAU 3.2 – Average, 5% percentile and standard-deviation of the Kou model parameters estimated from the option data between 2006 and 2011</caption><thead><tr><th></th><th>Average</th><th>5% percentile</th><th>Std-Dev</th></tr></thead><tbody><tr><td>p</td><td>0.64</td><td>0.84</td><td>0.24</td></tr><tr><td>η<sup>+</sup></td><td>0.16</td><td>0.28</td><td>0.06</td></tr><tr><td>η<sup>-</sup></td><td>0.15</td><td>0.28</td><td>0.07</td></tr><tr><td>λ</td><td>0.62</td><td>2.44</td><td>0.12</td></tr><tr><td>σ</td><td>18.29%</td><td>29.64%</td><td>0.08</td></tr></tbody></table>

tions. An example of the result on the calibration is shown in Figure 3.4.

A few remarks on the calibration can be made:

- Since the upward-sloping part of the smile is very small, the positive jumps are hardly calibrated in a reliable manner. However, the pricing of the gap option (section II.2.2) only needs the negative jumps intensity (i.e the downward-sloping part of the smile).

- The calibration is better on close-to-maturity options (as mentioned in [152]). It allows a better capture of instantaneous jump.

- The calibrated parameters will be used for hedging gap risk in the last section.

Figure 3.5 compares different discrete rebalancing frequencies with a jump modeling:

- Even for daily rebalancing, breaching the floor is unavoidable with the same probability as the two other frequencies.

- The three rebalancing frequencies give similar results when taking transaction costs into account.
<table><caption>TABLEAU 3.3 – Performance of the CPPI for different models and rebalancing frequencies</caption><thead><tr><th rowspan="2"></th><th colspan="3">Kou model</th></tr><tr><th>Daily</th><th>Weekly</th><th>Monthly</th></tr></thead><tbody><tr><th>Mean</th><td>146.28</td><td>147.10</td><td>147.57</td></tr><tr><th>Std-Dev</th><td>52.84</td><td>52.93</td><td>53.11</td></tr><tr><th>95% quantile</th><td>92.19</td><td>92.21</td><td>92.03</td></tr><tr><th>99.5% quantile</th><td>59.38</td><td>59.08</td><td>59.23</td></tr><tr><th>5% quantile</th><td>238.13</td><td>238.67</td><td>239.41</td></tr><tr><th>0.5% quantile</th><td>349.41</td><td>350.92</td><td>350.37</td></tr><tr><th>Rebalancing cost</th><td>0.92</td><td>0.45</td><td>0.26</td></tr></tbody></table>

The previous illustrations show that both the frequency of the rebalancing and the modeling affect the final value. The two metrics previously defined for different modeling assumptions

— The local probability of breaching the floor :

$$P_{t_i, t_{i+1}}^{\text{LBF}} := \mathbb{P} (V_{t_{i+1}} \le F_{t_{i+1}} | V_{t_i} > F_{t_i}).$$

— The overall probability of breaching the floor :

$$P^{\text{BF}} := \mathbb{P}(\exists t \in [0, T] : V_t \le F_t) = \mathbb{P}(V_T \le F_T).$$

— For Black-Scholes model in discrete-time rebalancing :

$$P_{t_i, t_{i+1}}^{\text{LBF}} = \mathcal{N} \left( - \frac{\log\left(\frac{m}{m-1}\right) + \left(\mu - r\right) \frac{T}{N} - \frac{1}{2} \sigma^2 \frac{T}{N}}{\sigma \sqrt{\frac{T}{N}}} \right).$$

and

$$P^{\text{BF}} = 1 - \prod_{i=0}^{N-1} (1 - P_{t_i, t_{i+1}}^{\text{LBF}}).$$

— For Kou jump process in continuous time :

$$P^{\text{BF}} = 1 - \exp \left( -T p \lambda \left( 1 - \frac{1}{m} \right)^{\frac{1}{\eta}} \right).$$

— Results depend on the model parameters and discretization time step :

— Gap risk goes to 0 as the rebalancing tends to be more frequent

— When considering a discontinuous path (jump models), even in continuous rebalancing the gap risk value > 0

## Impact of the ratchet feature

The ratchet feature is used by insurance companies to attract investors as it periodically locks in profit, see [43] and [9] for more details : at anniversary dates the guarantee is set to the highest value so

<footnote n="1">1. Implied volatilities are collected from Bloomberg dataset. The calibration is performed weekly.</footnote>
<table><caption>TABLEAU 3.4 – The probability of breaching the floor for different models and rebalancing frequencies.</caption><thead><tr><th>Model</th><th>Frequency</th><th>P<sup>BF</sup></th></tr></thead><tbody><tr><td rowspan="3">B&S</td><td>Monthly</td><td>9,07 × 10<sup>-5</sup></td></tr><tr><td>Weekly</td><td>1,2 × 10<sup>-10</sup></td></tr><tr><td>Daily</td><td>~ 0</td></tr><tr><td>Kou</td><td>Continuous</td><td>0,000410</td></tr></tbody></table>

far. The guarantee G becomes a time dependent function.

$$G_t = \begin{cases} V_0 & \text{if } t = 0 \\ \max(G_{t_{k-1}^*}, V_{t_k^*}) & \text{if } t = t_k^* \\ G_{t_k^*} & \text{if } t \in (t_k^*, t_{k+1}^*). \end{cases}$$

The bond floor is then defined as $F_t = G_t e^{-r(T-t)}$.

This feature has advantages and drawbacks. Locking-in the cash will ensure a higher guarantee but also reduces the cushion, the risky exposure and thus the upside potential risk.

The main results from figure 3.6 are :

- The mean and standard deviation of the final value increase with the rebalancing frequency (see table 3.5). This is justified by the path dependency of the guarantee which has a larger distribution with higher rebalancing frequency.

- The quantiles on the two tails of the final value distribution increase with the rebalancing frequency, while the distribution is shifted to the right with narrower body.

<table><caption>TABLEAU 3.5 – Final value metrics : Comparison between a CPPI without and with the ratchet feature</caption><thead><tr><th rowspan="2"></th><th colspan="3">Without ratchet</th><th colspan="3">With ratchet</th></tr><tr><th>Daily</th><th>Weekly</th><th>Monthly</th><th>Daily</th><th>Weekly</th><th>Monthly</th></tr></thead><tbody><tr><th>Mean</th><td>123.82</td><td>124.26</td><td>124.17</td><td>145.46</td><td>143.01</td><td>134.03</td></tr><tr><th>Std-Dev</th><td>41.96</td><td>43.29</td><td>47.25</td><td>100.08</td><td>81.75</td><td>45.60</td></tr><tr><th>95% quantile</th><td>99.99</td><td>99.68</td><td>90.88</td><td>100.61</td><td>100.53</td><td>99.99</td></tr><tr><th>99.5% quantile</th><td>99.99</td><td>97.57</td><td>77.84</td><td>99.99</td><td>99.94</td><td>98.27</td></tr><tr><th>5% quantile</th><td>214.18</td><td>216.42</td><td>222.23</td><td>268.74</td><td>261.73</td><td>219.52</td></tr><tr><th>0.5% quantile</th><td>289.15</td><td>292.33</td><td>314.58</td><td>700.97</td><td>603.28</td><td>359.59</td></tr><tr><th>P<sup>BF</sup></th><td>0.11</td><td>0.47</td><td>0.64</td><td>0.11</td><td>0.48</td><td>0.84</td></tr></tbody></table>

Consider the stopping time $\tau$ as the first time the portfolio value breaches the floor which does not depend on the bond floor level. The distribution of $\tau$ is the same in case of adding the ratchet, i.e. the probability of breaching the floor is not usually affected by the ratchet feature in theory. However, in our simulations, this probability in higher for the monthly rebalancing. This might be

## 3.3 Mitigating downside risk : Preventing from breaching the floor

### 3.3.1 Adjusting the multiplier to market conditions

By focusing on managing returns in downside markets, CPPI effectively manages portfolio volatility. Over the 5-year data (which included one bullish market, one bear market and a recovery), the CPPI strategy resulted in a slightly lower return – but also a significantly lower volatility. Additionally, the
worst one year return for the CPPI strategy was significantly less than that of the Index Portfolio.

The manager usually sets the multiplier at the beginning of the period. The risky exposure depends then on the evolving cushion. As the probability of breaching the floor may surge in market crash, or the manager might miss the subsequent market recovery, the multiplier needs to be adjusted accordingly with the market conditions.

A first approach to define a dynamic multiplier is the choice of the optimal m, deduced from the closed form solutions for optimal payoffs, and optimal certainty equivalent returns (CERs) using HARA utilities and log-normal distribution, see [142]. The authors give the following formula $m^* = \eta(\mu - r)/\sigma^2$ ($\eta$ here is the investor's sensitivity of risk tolerance to wealth). A particular case is the growth optimal leverage with $\eta = 1$ which is resulted in optimizing the growth rate of the leveraged strategy (cushion).

An alternative to the optimal multiplier is a Value-At-Risk based multiplier where investors choose the confidence level according to their risk tolerance as well but focused on tail risks.

Based on the weight $w_t^R$ of the Value-At-Risk Based Portfolio Insurance (VBPI) introduced by [106], and the expression of the risky exposure in both strategies ($E_t = m_t C_t = w_t^RV_t$), the expression for the multiplier at time t is:

$$m_t = \frac{1}{1 - \exp((\mu - r - \frac{1}{2}\sigma^2)(T - t) - z_p\sigma\sqrt{T - t})}$$

Since the dynamic multiplier depends on both volatility and return estimates, in order to improve its efficiency, $\mu$ and $\sigma$ can be re-estimated at each time step. However, the estimation of the drift is hardly accurate. Therefore, we will restraint the time dependency to the volatility. It will be re-estimated through a 3-month sliding window to take into account different market regimes.

The two approaches offer an interesting alternative to constant multiplier which lacks flexibility to market conditions. The comparison between these two approaches through a backtesting from 2006 to 2011, is illustrated in Figure 3.7. The focus on two periods (2006-2007 and post 2008 crisis), in Figure 3.8, illustrates that the VaR-based multiplier can perform better than the "optimal" one in bullish market and recovery (e.g 18% return Q2-2009 until Q1-2011 vs 11% in the post 2008 crisis). In contrast, during bear market, using the "optimal" multiplier (through $m < 1$) helps keep a relatively higher cushion but misses the recovery as it doesn't allow a high leverage.

In order to allow to participate in the market recovery to a greater extent, the multiplier is adjusted with a modified volatility estimator, either through a short-term exponentially weighted moving average (EWMA with $\lambda = 0.94$) realized volatility or an estimator based on implied volatility (of the strike consistent with the latest market returns). For example, if the underlying jumped 5% downward, the implied volatility with strike 95% will be chosen. For unavailable strikes, we use a linear interpolation. This strategy starts reinvesting into the risky asset as soon as Q3 2009, resulting in a higher performance by allowing the portfolio to capture more of the upside return when markets rebound. The backtesting in Figure 3.10 illustrates that the new multiplier is more reactive when adjusting with the implied volatility estimator. However, the 3-month realized volatility provides a higher multiplier and, when considering transaction costs, leads a lower cost of management.

Finally, the fixed frequency rebalancing is switched to a trigger rebalancing which occurs when the multiplier is out of a specific range chosen by the portfolio manager. In our case, on average the rebalancing frequency becomes every other day, which is consistent with the usual practice in CPPI asset management. At the same time, the cost of rebalancing is cut by half in comparison to a daily rebalancing (i.e. as low as a weekly or monthly rebalancing). Figure 3.11 illustrates the increasing performance
specifically under a range-bound high volatility regime, e.g. Q1-Q3 2008.

Adjusting the multiplier dynamically allows it to be more reactive to market conditions and explicitly dependent on the investor's risk aversion. However, it does not totally annihilate the downside risk in case of sudden jumps, where options may be useful to hedge those gap risks.

### 3.3.2 Hedging gap risks

The CPPI methodology will not necessarily protect the portfolio against a “black swan” event (such as a market crash of 20% in one day). To the extent that asset allocation shifts are implemented via underlying funds, the rebalancing trade can only occur at the end-of-day NAV. Even if futures are used to implement shifts intra-day, there can be gap movements in the futures markets. This is where a small gap risk protection sleeve can add value to the portfolio. To protect against such a “black swan” event, it is important to already have put options on market indices in the portfolio.

#### Vanilla Put option

A simple hedging strategy for the CPPI through embedded option can be constructed using short maturity put options. Touching the bond floor is mathematically equivalent to the cushion becoming negative. Assuming the event hasn't occurred up to time $t_k$, using equation (3.1), we have:

$$C_{t_{k+1}} < 0 \Leftrightarrow m \frac{S_{t_{k+1}}}{S_{t_k}} - (m-1)e^{r_{\frac{T}{N}}^T} < 0.$$

Hedging this risk is equivalent to forcing this quantity to be positive. This can be done by buying a put option at each of the CPPI rebalancing period with strike $(1 - \frac{1}{m})e^{r_{\frac{T}{N}}^T} S_{t_k}$ and as a maturity the CPPI rebalancing frequency. To hedge the whole portfolio the manager needs a number of $m \frac{C_{t_k}}{S_{t_k}}$ puts, which is the risky asset exposure. The discounted payoff in this case is $e^{-r_{\frac{T}{N}}^T} C_{t_k} ((m-1)e^{r_{\frac{T}{N}}^T} - m \frac{S_{t_{k+1}}}{S_{t_k}})^+$. The hedging cost at time $t_k$ can be written as:

$$\text{Cost}_{t_k} = m \frac{C_{t_k}}{S_{t_k}} \mathbb{E} \left[ \left[ \left(1 - \frac{1}{m}\right) e^{r_{\frac{T}{N}}^T} S_{t_k} - S_{t_{k+1}} \right]^+ \right].$$

Two approaches can be considered:

- The hedging costs (put prices) are deducted only afterwards from the portfolio value (which allows an estimation of how much the hedge would cost). In this case, the cushion follows the recursive relation:

$$C_{t_{k+1}} = C_{t_k} \left( m \frac{S_{t_{k+1}}}{S_{t_k}} + (1-m) e^{r_{\frac{T}{N}}^T} \right)^+$$

The cost of hedging can be computed as the sum of all put options prices necessary for the hedging:

$$C = \sum_{k=0}^{n-1} m \frac{C_{t_k}}{S_{t_k}} \mathbb{E} \left[ \left[ \left(1 - \frac{1}{m}\right) e^{r_{\frac{T}{N}}^T} S_{t_k} - S_{t_{k+1}} \right]^+ \right].$$

- In practice, the price of the puts used for the hedge will be deducted from the portfolio value at each step. This is translated in the second approach where the cushion dynamics follows the
recursive equation :

$$
\tilde{C}_{t_{k+1}} = e^{-r_{\frac{T}{N}}^T} \tilde{C}_{t_k} \frac{\left( m \frac{S_{t_{k+1}}}{S_{t_k}} + (1-m)e^{r_{\frac{T}{N}}^T} \right)^+}{\mathbb{E}^{\mathcal{Q}} \left[ \left( m \frac{S_{t_{k+1}}}{S_{t_k}} + (1-m)e^{r_{\frac{T}{N}}^T} \right)^+ \middle| \mathcal{F}_{t_k} \right]}.
$$

In Figure 3.13, we compare the effects of hedging using puts. We can notice the following :

- The guarantee is ensured and the manager no longer holds the risk of breaching the floor. However, once the put is exercised and the floor recovered, the manager needs to monetize in order to keep the guarantee until maturity.

- In terms of distributions, the CPPI distribution with a put hedging is a truncation of the classical CPPI where losses are cut (left tail limited by the guarantee).

Gap put option

An alternative risk mitigating action lies in the use of gap options which allow for a protection against sudden significant and persistent downside market moves : if a gap event occurs between two consecutive dates, the buyer receives the difference between the performance of the risky asset at gap $r = \frac{S_t}{S_{t-1}} - 1$ and the threshold J. In case of the CPPI, the proposed solution is a gap put option whose notional is the risky exposure with strike $J = 1/m$, where $m$ is the multiplier. We give in the following the main results for pricing of a gap option from [152]. We let the reader refer to the main article for details.

Suppose that the time to maturity T of a gap option is subdivided onto N periods of length h (e.g. days) : $h = \frac{T}{N}$. The return of the kth period will be denoted by $R_k^\Delta = S_{kh}/S_{(k-1)h}$.

Let $\alpha$ denote the return level which triggers the gap event and $k^*$ be the time of first gap expressed in the units of $h : k^* := \inf\{k : R_k^h \le \alpha\}$. The gap option is an option which pays to its holder the amount $f(R_k^{h*})$ at time $hk^*$, if $k^* \le N$ and nothing otherwise.

Assuming a deterministic interest rate $r$ and an i.i.d log returns $(R_k^h)_{k=1}^N$ and denote the distribution of $\log(R_1^h)$ by $p_h(dx)$. Then the price of a gap option is given by :

$$
G_h = e^{rh} \int_{-\infty}^{\beta} f(e^x) p_h(dx) \frac{1 - e^{-rT} (\int_{\beta}^{\infty} p_h(dx))^N}{1 - e^{-rT} \int_{\beta}^{\infty} p_h(dx)},
$$

with $\beta := \log(\alpha) < 0$.

It is complicated to obtain numerical results using this expression. Therefore, an approximate formula is used.

Let us assume $S_t = S_0 e^{X_t}$, where X is a Lévy process. Considering the hypothesis $rh \sim 10^{-4}$ and $h \to 0$, the following formula is obtained :

$$
G_h \approx \int_{-\infty}^{\beta} f(e^x) \frac{1 - e^{-rT} \left( \int_{\beta}^{\infty} p_h(dx) \right)^N}{1 - e^{-rT} \int_{-\infty}^{\beta} v(dx)}
$$

Assuming a Kou model (for its tractability and simplicity in integration) and considering the put payoff (i.e $f(x) = (K - x)^+$). The price then becomes :

$$
G_h \approx \frac{\lambda p \eta^{-}}{1 + \eta^{-}} K^{1+1/\eta^{-}} \frac{1 - e^{-T(r+\lambda p e^{\beta/\eta^{-}})}}{r + \lambda p e^{\beta/\eta^{-}}}
$$
with $p$ the probability that a given jump is negative, $\eta^-$ its intensity and $\lambda$ the poisson process rate. Moreover, for the CPPI we are interested in the payoff $((m-1)e^{rh} - m \frac{S_{kh}}{S_{(k-1)h}})^+$ which is equivalent to $(\frac{m-1}{m}e^{rh} - x)^+$ and thus, $K = (1 - 1/m)e^{-rh}$

The gap put option allows to cut the loss compensates for the loss as the portfolio value breaches the bond floor. However, insurance investors holding a CPPI who want to hedge it with gap option may face the following issues :

- The price of the gap option is usually sold higher than its theoretical cost for several reasons :

- The cost of the hedging the gap option for the bank may be quite higher because of the illiquidity of deep out of the money options that replicate it.

- The replicating formula is tricky to implement and interpret, as significantly model dependent (jumps multiple parameters, lack of robustness).

- Actually, the gap option proposed by the bank might have a different design and payoff from the one considered for the hedge.

- The bank usually hedges the gap up to the first order only.

- The gap risk is borne by the bank only if there is some reconciliation by the insurance company within 24/48 hours, out of which the insurer bears herself the gap risk. As a result, operational risks are significant and represent a major part of the economic capital requirements (e.g. under Solvency II framework).

<table><caption>TABLEAU 3.6 – Final value metrics : Comparison between different hedging strategies</caption><thead><tr><th rowspan="2"></th><th colspan="3">Hedging strategies</th></tr><tr><th>Vanilla Put Hedge 1</th><th>Vanilla Put Hedge 2</th><th>Gap Option</th></tr></thead><tbody><tr><th>Mean</th><td>136.97</td><td>133.35</td><td>134.98</td></tr><tr><th>5% quantile</th><td>218.70</td><td>215.40</td><td>217.00</td></tr><tr><th>0.5% quantilee</th><td>277.21</td><td>273.53</td><td>275.22</td></tr><tr><th>Hedging cost</th><td>N.C</td><td>2.26</td><td>1.08</td></tr></tbody></table>

## 3.4 Conclusion

In this article we have presented a study of the CPPI as an insurance contract, a review of its theory and practice as well as its modeling and hedging issues for a risk/return/cost perspective. The main conclusions are :

- Continuous CPPI is only theoretical : given market frictions and the probability of not ensuring the guarantee, all the more that jumps occur more than not.

- As a result, jump processes are a valuable input for the CPPI modeling : they allow to catch a probability of breaching the floor different than zero (even in the continuous-time framework; [83] and [84] came up with the same conclusion) and therefor, detect, define and hedge gap risk.

- Correctly choosing and adjusting the multiplier dynamically significantly reduce the downside risk according to a Value-At-Risk indicator : The multiplier decreases in period of turmoils reducing the risky exposure and increases back during market recovery.

- Hedging the gap risk is possible through two types of options : Vanilla Puts and Gap Put options. The first one is more common due to liquid assets, but the hedging cost may turn out to be too expensive and the maturity too limited. The second type of options is less liquid (bought only through an agreement) but is cheap.
<figure><img src="image_32.png" /><figcaption>FIGURE 3.1 – Evolution of an investment in the S&P500 for the period Q1 2006 to Q4 2010</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>FIGURE 3.2 – Backtesting and distribution of the three various rebalancing frequencies under B&S model.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>FIGURE 3.3 – Comparison between the Buy & Hold strategy, CPPI with m = 3 and CPPI with m = 6 through backtesting (S&P500)</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>FIGURE 3.4 – Calibration of the Kou model using 1-month maturity call options price on the S&P500</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>FIGURE 3.5 – Simulation and distribution of the three various rebalancing frequencies under Kou model.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>FIGURE 3.6 – The figure on the top is a backtesting on the previous set of data to compare a classical iCPPI and one with the ratchet feature. The three histograms on the bottom are those of the final value distribution for the three different rebalancing frequencies on the iCPPI with ratchet.</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>FIGURE 3.7 – Comparison between different multipliers (VaR-based with $p = 99.5\%$ and the optimal one with risk tolerance $\eta = 0.2, 0.4$ and $1$) based on Realized Volatility</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>FIGURE 3.8 – Focus on two bullish market periods where the CPPI with the VaR-based m performs better than the optimal one</figcaption></figure>
<figure><img src="image_40.png" /><figcaption>FIGURE 3.9 – Comparison between dynamic multiplier based on RV and on IV through backtesting</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>FIGURE 3.10 – Comparison between dynamic multiplier based on RV and on EWMA through backtesting</figcaption></figure>
<figure><img src="image_42.png" /><figcaption>FIGURE 3.11 – Comparison between trigger rebalancing vs fixed frequency rebalancing</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>FIGURE 3.12 – Comparison between the dynamic multiplier and an adjusted one based on a manager decision depending on market recovery</figcaption></figure>
<figure><img src="image_44.png" /><figcaption>FIGURE 3.13 – Comparison between no hedging and put hedging in its two approaches</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>FIGURE 3.14 – Comparison between a vanilla and a gap option hedging</figcaption></figure>
Deuxième partie

**Market impact and volatility**
<figure><img src="image_46.png" /></figure>
# Chapitre 4

## On Optimal Options Book Execution
Strategies with Market Impact

**Abstract**— We consider the optimal execution of a book of options when market impact is a driver of the option price. We aim at minimizing the mean-variance risk criterion for a given impact function. First, we develop a framework to justify the choice of our impact function. Our model is inspired from Leland's option replication with transaction costs. The option effective price is then expressed through a Black-Scholes like PDE with a modified volatility that depends on the size of the trade. We set up a stochastic control framework and solve an Hamilton-Jacobi-Bellman equation using finite differences methods. The simple expected cost problem suggests that the strategy is characterized by a convex increasing trading speed, in contrast to the equity case where the optimal strategy results in a constant trading speed. However, in this framework, the underlying price does not seem to affect the agent's decision. By taking the agent risk aversion into account through the variance, the strategy seems to be more sensitive to the underlying price evolution, urging the agent to trade faster at the beginning of the strategy.

**Keywords** : Market impact; option pricing; optimal execution; mean-variance; stochastic control; HJB equation

### 4.1 Introduction

This chapter addresses the optimal execution of a large portfolio of options, from a starting composition to a specified final one, within a specified period of time. This problem was first introduced by Bertsimas and Lo [27] in the context of a large block of equity portfolio. They define best execution as the dynamic trading strategy that minimizes the expected cost of trading over a fixed period. Given a fixed block of shares and a price-impact function that drives the execution price of an individual trade – that depends on the traded shares and market conditions –, they obtain the optimal *sequence* of trades. Their approach focuses on the expected cost and ignores the volatility of revenues for different trading strategies. Later on, Almgren and Chriss [7] work on a more general framework, using the variance as a suitable penalty for the uncertainty of the cost. They developed the fundamentals of the optimal execution under *market impact* constraints. Seventeen years later, this problem is still attracting the interest of both practitioners and academics. Since then, there was a tremendous growth in optimal trading algorithms, fueled by high frequency trading and powerful computer units which increased market efficiency and decreased the free-lunch opportunities.
4.1.1 Motivation

Derivatives can be used to hedge or to speculate, and many financial institutions include them as a considerable part of their books, occupying a large proportion in their trading activity. Options are common derivatives that allow making profit when market goes down or even sideways. However, the most useful function of options remains the hedging against undesirable movements of an asset. Let us take the example of an insurance company. Insurance liabilities are characterized by three main features : long term time horizon, large volumes and significant market exposure. Given the market uncertainties, the insurance company needs protection strategies in order to hedge its significant risk exposure. In that respect, the use of put options is very common. Buying such hedging portfolios on a significant scale requires taking into account the size of the transaction explicitly.

There is a large literature which develops various parametric models of derivative prices according to the no-arbitrage theory, see [31] and [126]. Bates who surveyed this literature in [17], emphasizes that it cannot fully capture the empirical properties of option prices. He concludes that there is a need for a new approach to pricing derivatives that focuses on the "financial intermediation of the under- lying risks by option market-makers", see [17]. In [34], Bondarenko addresses the "overpriced puts puzzle" through studying the historical prices of the S&P 500 put options. He finds that their price has been too high and incompatible with canonical asset-pricing models. Such mismatch is illustrative of a potential equilibrium premium stemming from supply/demand imbalances, which is not explicitly considered by traditional models. Trying to solve this dilemma, Gârleanu and co-authors, see [85], conduct an empirical work to develop a demand-based theory. In their approach, market makers who incur higher unhedgeable risks will move the price up if the net demand is positive and down if it is negative. Even though the equity options market is in zero net supply, negative net demand by the end users in the equity options market, can be viewed as economically equivalent to the positive net supply in the equity market. Another empirical study was conducted in [52] to determine the effect of option and stock illiquidity on delta-hedged equity option returns. They present strong evidence on illiquidity premium in option markets, using effective spreads on a large number of underlying firms with intraday trades and quotes.

The theoretical literature also studied the impact of illiquidity of the underlying stocks on option prices. In a frictionless, complete-market model, the price of the option can be replicated by trading in the underlying asset and risk free bond. If the asset is subject to additional transaction/market impact costs, this should affect the return on options. In [115], Leland provides a theoretical analysis of this effect using a hedging argument. Because option market makers are net long in equity option markets, they need to create a synthetic short option using the underlying stock. Another substantial amount of literature on option pricing and hedging in the presence of market impact has been developed. We mention, among others, [1, 118, 119] in which the problem of option hedging for a large trader who experiences market impact is considered. An Hamilton-Jacobi-Bellman (HJB) equation is derived and a fully non-linear pricing partial differential equation (PDE) is solved. This is known as the feedback effect and yields to an option pricing with the adjusted underlying price. Other references on the subject deserve to be mentioned, see [77, 78, 117].

To our best knowledge, the optimal execution of a portfolio of options under market impact constraints has been forsaken, and our work is the first one to address this problem. Our paper seeks to incorpo- rate the market impact on options into the framework of optimal execution. To do so, the Leland's transaction costs approach will be used to derive the option effective price, and the Almgren expected mean-variance optimization framework will be adopted to define what is the best execution strategy.
### 4.1.2 A review of the linear market impact model for equity

A market impact model aims at describing the quantitative feedback of a large order, one that affects the execution price and quantifies illiquidity.

In [5] and [7], the price impact is presented as a combination of two components:

- A permanent component that reflects the information transmitted to the market by the buy/sell imbalance.

- A temporary component that reflects the price concession needed to attract counter-parties within a specified short time interval.

We consider an asset S which is traded continuously. The number of shares of the traded asset is described by an absolutely continuous trajectory $t \to x_t$, $\dot{x}_t$, its derivative w.r.t time corresponds to the speed of trading of the security. In the absence of market impact, the asset is modeled by a geometric Brownian motion (GBM)

$$S_t = S_0 e^{-\frac{1}{2}\sigma^2 t + \sigma W_t} \quad (4.1)$$

When market impact is taken into account, the execution price $\tilde{S}_t$ is defined by

$$\tilde{S}_t = S_t(1 + \eta\dot{x}_t + \gamma(x_t - x_0)), \quad (4.2)$$

where S is the unaffected stock price process, and $\eta$ and $\gamma$ are constants.

The term $\eta\dot{x}_t$ corresponds to the temporary or instantaneous impact of trading $\dot{x}_t dt$ shares at time $t$ and only affects this current order. The term $\gamma(x_t - x_0)$ is the permanent price impact which was accumulated by all transactions until time $t$.

Both impacts in the previous model are linear on the trading volume. In practice, however, even though strong evidence from [100] argues that permanent impact is linear, the linearity of the instantaneous impact is an unrealistic assumption. Perold and Salomon [141] argue that the liquidity premium per share demanded by the market will be either a convex or a concave function of the block size. Other independent empirical studies have demonstrated that the price change induced by the sequential execution of a total volume X follows an approximate $\sqrt{X}$ law, see [16, 44], or of a 3/5 power law, see [8]. In this paper, we focus exclusively on linear price-impact function since our main interest is the optimal execution strategy.

### 4.1.3 Main results and organization of the paper

In Section 4.2 we derive a put pricing model that takes into account market impact using Leland's transaction costs framework, see [115]. We consider the point of view of a market-maker who sells put options and immediately takes position on the underlying asset to hedge his options. The hedging comes with an additional cost due to trading the asset. The additional cost is then incorporated into the option price. Our option execution price can be expressed through a Black-Scholes like PDE with an enlarged volatility such that

$$\tilde{\sigma}^2 = \sigma^2 + f(t, \dot{x}_t, x_t, \sigma),$$

where $\sigma$ is the asset volatility and $f$ is the market impact function (depends on time, volatility, inventory and trading speed).

We can either consider the PDE directly or express the put price using the Black-Scholes closed formula with the enlarged volatility when $\sigma$ is constant. We can also extend the formula using a Taylor
approximation w.r.t the volatility. By doing so, we can recover an additive formulation of the impact function. Either way, we conclude that the impact on the option drives its volatility and that the impact is a function of the option vega.

Once given the put price in presence of market impact, we will set up the optimal execution framework in Section 4.3. We consider an agent who has a fundamental need for option exposure. Such agent is denoted as "end-user" and wishes to acquire a large quantity of put options but needs to minimize the cost of its acquisition. The solution to such optimization leads to a closed formula in the expected cost minimization framework. The trading speed in this case is of the inverse of a 3/2 power law, in contrast to the equity case which is known to result in a constant trading rate strategy.

In Section 4.4 we choose the mean-variance framework as in previous works, see [6, 7, 121]. An approximation to the variance is performed as its drift part is complicated to compute explicitly. Such approximation was already used in the Almgren-Chriss framework. Taking the impact as an additional term in the option price process, we find an HJB equation which leads to a quasi-linear PDE. Using a proper parametrization of the state variable, we are able to reduce the dimension of the problem and solve the PDE numerically through finite differences methods in Section 4.5. Results are presented in Section 4.5.2

Finally, in Section 4.6 we use a semi-Lagrangian approach for the mean-variance framework when the option price follows the PDE formulation. References on this method can be found in [49, 64, 75]. This allows to use a local volatility model and consider permanent impact as well.

## 4.2 Market impact model: A transaction costs approach

The main objective of the following section is to express the option execution price in the presence of market impact. We build our model upon two areas of the cited literature. Market impact in equity market, and option pricing under transaction costs. The main result of the section is the option price PDE expressed through Formula (4.5).

### 4.2.1 The transaction costs approach revisited

We consider a market with two agents; the "option market-maker" (or sell-side trader) who is issuing the vanilla put option and dynamically hedging his position, and the "end-user" (or buy-side trader) who buys the option in order to hedge her fundamental risk. If end-users were able to hedge perfectly – as in a Black-Scholes framework – then option prices are determined by the classical no-arbitrage theory without market impact constraints. In reality, however, perfect hedging is not possible (because of the impossibility of trading continuously, stochastic volatility, jumps, transaction costs, etc...).

In his seminal paper [115], Leland suggests a modified Black-Scholes approach to a contingent claim pricing with proportional transaction costs. Fortunately, we can assimilate market impact to transaction costs. We exploit this framework in order to develop an option pricing model that allows to incorporate the traded quantity as an additional variable.

The dynamic hedging of the option consists in investing a proportion $\Delta$ on the underlying and the remaining on a non-risky asset. Let $\hat{T}$ be the maturity of the option, issued at time 0, and let us consider a discrete grid with $n$ revision dates. The non-risky asset is the *numeraire* $S_0$ and the risky asset is
given under the martingale measure by the SDE

$$dS_t = \sigma S_t dW_t, \quad 0 \le t \le \hat{T}.$$

From the option market-maker perspective, the trading involves proportional impact rate $I_0$, fixed at the date on which she sells the option (and starts delta-hedging it). This is similar to a transaction cost in the Leland framework. The current value of the portfolio process at time $t$ is defined by

$$V_t^n = V_0^n + \int_0^t \Delta_u^n dS_u - \sum_{t_i \le t} I_0 S_{t_i} |\Delta_{i+1}^n - \Delta_i^n|, \quad t \le \hat{T} \quad (4.3)$$

where:

- $t_i = t_i^n = i/n, 0 \le i \le n, t_0 = 0, t_n = \hat{T}$ are the revision dates.

- $\Delta^n = \Delta_i^n$ on the interval $]t_{i-1}, t_i]$, $\Delta_{n+1}^n := \Delta_i^n$.

- $\Delta_i^n$ is $\mathcal{F}_{t_{i-1}}$-measurable.

$\Delta^n$ corresponds to the trading strategy. The number of shares of the risky asset that the holder possesses in the period $i$ is then $\Delta_i^n$. The dynamics (4.3) means that the portfolio process $V^n$ is self-financed and in presence of market impact (which is proportional to the traded volume).

When market is complete without friction, the option price is exactly replicated by the terminal value of the self-financing portfolio:

$$V_t = \mathbb{E}[(K - S_{\hat{T}})^+ | \mathcal{F}_t] = P(t, S_t) = V_0 + \int_0^t \partial_S P(u, S_u) dS_u, \quad t \le \hat{T},$$

where $P$ is solution of the PDE:

$$\begin{cases} \partial_t P(t, S) + \frac{1}{2} \sigma^2 S^2 \partial_{SS} P(t, S) = 0, & (t, S) \in [0, \hat{T}[\times]0, \infty] \\ P(\hat{T}, s) = (K - s)^+, & s \in ]0, \infty[. \end{cases} \quad (4.4)$$

Following Leland's approach, we construct a strategy which can be perceived as a modified-Delta Black-Scholes replication formula. The constant volatility $\sigma$ is replaced by $\tilde{\sigma}$ in order to compensate for the market impact cost. The "enlarged volatility" $\tilde{\sigma}$ is defined by

$$\tilde{\sigma}^2 = \sigma^2 + \sigma I_0 n^{1/2} \sqrt{\frac{8}{\pi}}.$$

Therefore, the modified option price under market impact follows the PDE:

$$\begin{cases} \partial_t \tilde{P}(t, S) + \frac{1}{2} \tilde{\sigma}^2 S^2 \partial_{SS} \tilde{P}(t, S) = 0, & (t, S) \in [0, \hat{T}[\times]0, \infty] \\ \tilde{P}(\hat{T}, s) = (K - s)^+, & s \in ]0, \infty[. \end{cases} \quad (4.5)$$

Let us give, in the following, the intuition behind the result inspired from Leland strategy. Assuming the solution $\tilde{P}$ to (4.5) is smooth enough, we have

$$\tilde{P}(t, S_t) = \tilde{P}(0, S_0) + \int_0^t \partial_S \tilde{P}(u, S_u) dS_u + \frac{1}{2} \int_0^t [\sigma^2 - \tilde{\sigma}^2] S_u^2 \partial_{SS} \tilde{P}(u, S_u) du.$$

Therefore, $\tilde{P}(t, S_t)$ can be seen as the continuous version of a portfolio process (4.3), provided that $\Delta_i^n = \partial_S \tilde{P}(t_{i-1}, S_{t_{i-1}})$ and the drift term in the formula above corresponds to the cumulative market impact cost. We want to equate the two following increments:

$$\frac{1}{2} [\sigma^2 - \tilde{\sigma}^2] S_u^2 \partial_{SS} \tilde{P}(u, S_u) du,$$
and

$$-I_0 | \partial_S \tilde{P}(u + \Delta u, S_{u+\Delta u}) - \partial_S \tilde{P}(u, S_u) | S_{u+\Delta u}.$$

To do so, we use the Taylor expansion

$$\begin{aligned} \partial_S \tilde{P}(u + \Delta u, S_{u+\Delta u}) - \partial_S \tilde{P}(u, S_u) &\approx \partial_{S,t} \tilde{P}(u, S_u) \Delta u + \partial_{SS} \tilde{P}(u, S_u) (S_{u+\Delta u} - S_u) \\ &\approx \partial_{SS} \tilde{P}(u, S_u) (S_{u+\Delta u} - S_u), \end{aligned}$$

where

$$S_{u+\Delta u} - S_u \approx \sigma S_u (W_{u+\Delta u} - W_u).$$

Since $\partial_{SS} \tilde{P} \ge 0$, we should look for $\tilde{\sigma}$ such that

$$\frac{1}{2}[\sigma^2 - \tilde{\sigma}^2]\Delta u \approx -I_0 \sigma |W_{u+\Delta u} - W_u| \frac{S_{u+\Delta u}}{S_u}.$$

Then, considering the conditional expectation knowing $\mathcal{F}_u$, and the equalities

$$\mathbb{E} | W_{u+\Delta u} - W_u | = \sqrt{\Delta u} \sqrt{\frac{2}{\pi}}, \quad \frac{S_{u+\Delta u}}{S_u} = 1 + \sigma (W_{u+\Delta u} - W_u),$$

we deduce that

$$\tilde{\sigma}^2 = \sigma^2 + I_0 n^{1/2} \sqrt{\frac{8}{\pi}} \sigma. \qquad (4.6)$$

**Remarks 1. :**

- The mathematical justification of the approximate replication principle in the Leland transaction costs framework brought a lot of attention and turned out to be quite difficult to obtain. The replication principle fails to be true for constant transaction costs. In particular, the agent should choose and fix the number of revision dates $n$ large enough so that the mean square hedging error is controlled as there is a limit error as $n \to \infty$. In [122], Lott gives the first rigorous result : the approximation error tends to zero in probability if the transaction costs coefficient depends on $n$ and decreases on a rate proportion to $n^{-1/2}$ (in this case, $\tilde{\sigma}$ does not depend on $n$).

- To solve the previous issues, we will fix the number of revision dates to be of a reasonable order to provide a reasonable price. We do not aim at reviewing the foundations of the replication issues of the Leland framework, instead, our main interest is finding a simple and intuitive pricing formula for the option market-maker, one that takes into account additional fees due to market impact. As a consequence, the end-user will pay the option at a higher price than what the Black-Scholes pricing model suggests. The quantity of his trades will be an important factor.

- In the Black-Scholes framework, the implied volatility is constant, and so is the "enlarged volatility". In practice, however, the implied volatility is not constant, i.e depends on the maturity and the moneyness. One of the common market practices is to take local volatility models. Fortunately, the Leland approach is still valid in this case as showed in [117]. One just needs to replace the constant volatility $\sigma$ by a local volatility one $\sigma(t, S_t)$. The justification above remains the same.
4.2.2 The option price under market impact and the market impact function

In the following, we properly link the market impact framework to the transaction costs one. We consider a universe with the two agents; option market-makers and end-users. The end-user acquires, at time $t < \hat{T}$ (where $\hat{T}$ is the maturity of the option at time 0, and $\tau = \hat{T} - t$ is the residual time to maturity at time $t$), $x_t$ options with a trading speed $\dot{x}_t$. The option market-maker, on the other hand, meets the end-user's demand by selling the exact quantity at the exact same trading speed on the "option market". He then hedges these options on the "equity market" <fnref n="1" />. In doing so, he expects the asset price to be affected. The market impact (on the asset) is decomposed, as exposed in the introduction through Equation (4.2), into a temporary $I^{\text{temp}}$ and a permanent $I^{\text{perm}}$ impact, and both are linear on the quantity. At the infinitesimal time between $t$ and $t+dt$, the temporary impact resulting in selling $dx_t = \dot{x}_t dt$ options is expressed by

$$I_t^{\text{temp}} dt = \eta \dot{x}_t dt.$$

The permanent impact is a result of all previous trades up to time $t$. Note that at time $t$, the market-maker has already sold $(x_t - x_0)$ options and took position on the underlying asset to hedge them. Therefore, between $t$ and $t+dt$ we have

$$I_t^{\text{perm}} dt = \gamma(x_t - x_0) dt$$

The combination of the two impacts gives the total market impact rate at time $t$

$$I_t = \eta \dot{x}_t + \gamma(x_t - x_0). \quad (4.7)$$

Without adding complexity, the reasoning in Section 4.2.1 is straightforward in this case. We adapt the formula of the augmented volatility expression in Equation (4.6) at the arbitrary time $t \in [0, T]$, where $T < \hat{T}$ is the end time of the trading strategy. The revision step $h$ is fixed instead of the number of revision dates, regardless of the time-to-maturity $\tau = \hat{T} - t$ of the option. The number of revision dates becomes $n = \frac{\hat{T}-t}{h}$ and the augmented volatility at an arbitrary time $t$ can be expressed as

$$\tilde{\sigma}_t^2 = \sigma^2 + (\tilde{\eta}\dot{x}_t + \tilde{\gamma}(x_t - x_0))\sqrt{\hat{T} - t}\sigma, \quad (4.8)$$

where $\tilde{\eta} = \eta\sqrt{\frac{8}{h\pi}}$ and $\tilde{\gamma} = \gamma\sqrt{\frac{8}{h\pi}}$ are constants depending on the underlying impact factors $\eta$ and $\gamma$ and the revision step $h$ is fixed to be one day (meaning that the option market-maker will revise her hedging position on a daily basis). Finally, the option price under market impact $\tilde{P}$ follows the PDE:

$$\left\{ \begin{array}{l} \partial_t \tilde{P}(u, S) + \frac{1}{2} \tilde{\sigma}_t^2 S^2 \partial_{SS} \tilde{P}(u, S) = 0, \quad (u, S) \in [t, \hat{T}[\times]0, \infty] \\ \tilde{P}(\hat{T}, s) = (K - s)^+. \end{array} \right. \quad (4.9)$$

**Remark 6.** The term $\tilde{\gamma}(x_t - x_0)$ comes from previous trades (hedging options sold before 0 and $t$). By doing so, to the value of the "imaginary" asset price becomes $\tilde{S}_{t^-} = S_t + \gamma(x_t - x_0)S_t$. Then, by selling $\dot{x}_t dt$ at time $t$, and taking position of the underlying asset, the affected asset price becomes $\tilde{S}_t = \tilde{S}_{t^-} + \eta \dot{x}_t S_t$. Since the replication strategy at time $t$ is based on the asset price $S_t$ at time $t$, the impact term in the enlarged volatility is proportional to $\tilde{S}_t - S_t$, and therefore we have Equation (4.8).

**Proposition 1.** The execution price $\tilde{P}$ of an option sold at time $t$ can be written using Black-Scholes closed-formula with $\tilde{\sigma}_t$ as the volatility parameter:

$$\tilde{P}(t, S_t) = BS(t, S_t, K, r = 0, \hat{T}, \tilde{\sigma}_t),$$

<footnote n="1">1. We consider that the option market, in which the option is traded, is separated from the equity market, in which the option is hedged. This allows to justify the fact that the buy-side trader does not observe any changes in the asset price.</footnote>
where BS corresponds to the Black-Scholes closed-formula of the put option price, S_t is the option spot, K its strike, r the risk-free rate chosen to be 0, T the option maturity and σ the enlarged volatility defined by

$$ \tilde{\sigma}_t^2 := \sigma^2 + (\tilde{\eta}\dot{x}_t + \tilde{\gamma}(x_t - x_0))\sqrt{T-t\sigma} $$

Using a simple Taylor approximation to the first order around the spot variance V = σ², and assuming the impact term is negligible to σ², we can rewrite the above expression as a sum of the Black-Scholes option price with volatility σ and an additional term corresponding to the option market impact:

$$ \begin{aligned} \tilde{P}(t, S_t) &\approx P(t, S_t) + (\tilde{\sigma}_t^2 - \sigma^2) \partial_{\sigma} P(t, S_t) \\ &\approx P(t, S_t) + \frac{1}{2} \{ \tilde{\eta}\dot{x}_t + \tilde{\gamma}(x_t - x_0) \} \sqrt{T-t\sigma} (t, S_t), \end{aligned} \quad (4.10) $$

where v(t, S_t) = ∂_σ P is the Black-Scholes vega of the option, calculated on the constant volatility σ

$$ v(t, S_t) = \sqrt{T-t} S_t N'(d_1) = \sqrt{T-t} KN'(d_2), $$

where

$$ \begin{aligned} N'(x) &= \frac{1}{\sqrt{2\pi}} e^{-x^2/2} \\ d_1 &= \frac{\log \frac{S_t}{K} + \frac{1}{2}\sigma^2(\hat{T} - t)}{\sigma\sqrt{\hat{T} - t}} \\ d_2 &= \frac{\log \frac{S_t}{K} - \frac{1}{2}\sigma^2(\hat{T} - t)}{\sigma\sqrt{\hat{T} - t}} = d_1 - \sigma\sqrt{\hat{T} - t}. \end{aligned} $$

It follows that the impact function is defined by

$$ f(t, S_t, \dot{x}_t, x_t) := \tilde{P}(t, S_t) - P(t, S_t) = \frac{1}{2} \{ \tilde{\eta}\dot{x}_t + \tilde{\gamma}(x_t - x_0) \} \sqrt{T-t\sigma} (t, S_t). \quad (4.11) $$

**Remarks 2. :**

- The impact term on the underlying asset is not observed by the end-user. This assumption is made for the sake of deriving the market impact function on the option only. In fact, the sell-side trader bases her hedging portfolio on an "imaginary" asset whose price is subject to market impact. This effect will not be seen in the "real" asset and will not be perceived by the end-user. We can see it as if the market in which options are sold is separated from the market where the underlying asset is traded for the hedging. Such an assumption is not absurd given that the end-user is assumed to trade the option exclusively.

- Equation (4.10) allows to write the option execution price as the sum of martingale (the option price in the absence of market impact, i.e η = γ = 0) plus a positive term reflecting the additional cost due to market impact. The end-user trades the option exclusively and pays the additional hedging costs fixed by the sell-side trader and given by Equation (4.11). Of course the PDE (4.9) is more precise and allows to deal with both the Black-Scholes case and a local volatility model. However, the additive formula is more convenient to build an Almgren-like optimal execution framework.

- The vega-dependent impact function is a result of the agent adjusting the replicating portfolio according the "imaginary" price movement.
— The number of revision dates n can be fixed instead of the revision step h, which has the advantage of suppressing the term $\sqrt{\hat{T}-t}$ in the formula. However, by doing so, the sell-side trader needs to adjust the revision step h = $\frac{\hat{T}-t}{n}$ which decreases as time increases. The first option seems, however, more logical. The sell-side trader will revise her position on a daily basis for example (i.e h = 1/252).

— Using the Vega-Gamma relationship (i.e v = στS²Γ), we can rewrite the impact term as a function of the option Gamma:

$$f(t, S_t, \dot{x}_t, x_t) = \frac{1}{2} \{\tilde{\eta}\dot{x}_t + \tilde{\gamma}(x_t - x_0)\} \sigma S_t^2 (\hat{T} - t)^{3/2} \Gamma(t, S_t), \quad (4.12)$$

$$\text{where } \Gamma(t, S_t) = \frac{N'(d_1)}{S_t \sigma \sqrt{\hat{T}-t}} = \frac{KN'(d_2)}{S_t^2 \sigma \sqrt{\hat{T}-t}}.$$

**Proposition 2.** In the Black & Scholes framework, the put option effective price is written as the following:

$$\tilde{P}(t, S_t, \dot{x}_t, x_t) = P(t, S_t) + \frac{1}{2} \{\tilde{\eta}\dot{x}_t + \tilde{\gamma}(x_t - x_0)\} \sigma S_t^2 (\hat{T} - t)^{3/2} \Gamma(t, S_t), \quad (4.13)$$

where:

— $\tilde{\eta}$ controls the temporary impact strength in $ \times $ hour / N of options.

— $\tilde{\gamma}$ controls the permanent impact strength and is in $N$ shares.

— $x_t$ is the quantity held at time t and $\dot{x}_t$ is the speed of trading in number of options per time unit.

— $\Gamma$ is the delta sensitivity w.r.t to the asset price ($\Gamma(t, S_t) > 0$).

**Remark 7.** Buying the option (i.e $\dot{x}_t > 0$) will typically lead to increasing its price, hence the execution price $\tilde{P}_t \ge P_t$ for $t \in [0, T]$.

**Numerical experiment :**

We illustrate in Figure 4.1 the put price as a function of the moneyness for the Black-Scholes case. We seek to compare the price with and without market impact and justify the use of the Taylor expansion. We consider that the rate of trading is constant. The figure shows that the price when the option is subject to market impact is greater than the price without market impact. This is what we expect when we have a demand pressure on the option. On the other hand, the approximation by the Taylor expansion is almost equal to the closed formula. This will be very convenient to solve the optimal trade execution problem.
<figure><img src="image_47.png" /><figcaption>FIGURE 4.1 – Put price as a function of the moneyness.</figcaption></figure>

The Black-Scholes case results in a shift in the implied volatility (but a constant one in both cases). Solving PDE 4.9 with a Constant Elasticity Volatility model (CEV) gives the illustrations in Figure 4.2. In addition to a higher volatility across the moneyness, the increase in the implied volatility is larger for out of the money options.

<figure><img src="image_48.png" /><figcaption>FIGURE 4.2 – The CEV model volatility smile and put price evolution</figcaption></figure>
## 4.3 The optimal execution problem

In the previous section, we managed to define an effective price for an option under market impact. We assumed that the issuer prices the option by taking its replication price and an adjustment to his position when the underlying price is impacted by the order. We pointed out an important assumption; we described the impact on the asset as being "imaginary", meaning that the end-user does not see the asset price moving because of the strategy. The option price under market impact was found to verify a Black-Scholes like PDE with an "enlarged" volatility. The dependence to the position $x$ and trading speed $\dot{x}$ is linear and increasing (the higher the trading speed and quantity, the higher the volatility and thus the option price).

Particularly, an agent who is willing to trade a large quantity of options will see the impact as an important dilemma. If he trades rapidly, then the actual cost of her strategy will be different from the one she anticipated. In real life, the agent is exposed to price manipulation, liquidity issues and market impact, especially for large trades; the cost of placing one large order to close his position will be far greater than the sum of infinitely small orders differed in time.

Many works have been held, in the equity market, to show that because of the trades size of typical institutional investors, they are usually broken up into smaller ones and executed over the course of several days, see for example [48, 108]. Chan and Lakonishok [48] for example, show that only 20% of the market value of the trades splits in their set of data are completed within a day, and that over 53% are spread over four trading days or more. For this reason, best execution can not be defined as a single number in the context of a single trade. It is a strategy that unfolds over the course of several days and which ought to adapt to changing market conditions. This is intuitively true for market options as well, even though very little literature deals with this issue.

Taking this into account, instead of executing his orders at once, the agent has to split them over the time interval $[0, T]$ by means of a dynamic order execution strategy.

### 4.3.1 The general framework

Let us consider a trade execution strategy in which an initial long or short position of X options with fixed strike K and maturity $\hat{T}$ is liquidated by a fixed time horizon $[0, T]$, where $T < \hat{T}$ is the end time. We describe such a strategy by the asset position $x_t$ held at time $t \in [0, T]$. The initial position $x_0$ is positive for sell strategy and negative for buy strategy. The condition $x_{T^+} = 0$ assures that the initial position has been unwound by time T. The path $x = (x_t)_{t \in [0, T]}$ will be nonincreasing for pure sell strategy and nondecreasing for a pure buy strategy.

We restrict our framework to pure buy strategies. The end-user's purpose is to hedge the risk of a complex product (structured product, Variable Annuity, etc...) indexed on an underlying asset, by acquiring vanilla put options on that same underlying asset. Let $(\Omega, \mathcal{F}, \mathbb{P})$ be the usual probability space on the filtration $(\mathcal{F}_t)_{t \in \mathbb{R}_+}$ satisfying the usual assumptions. In the absence of market impact and under a null risk-free rate, the no-arbitrage price of a put option is defined by:

$$P_t = \mathbb{E}_{\mathbb{Q}} \left[ (K - S_{\hat{T}})^+ \mid \mathcal{F}_t \right]$$

where $\mathbb{E}_{\mathbb{Q}}$ is the expectation under the risk-neutral probability measure $\mathbb{Q}$ in which the asset price is a martingale. At each time $t$, $\dot{x}_t dt$ options are bought at price $\tilde{P}_t$. The effective price $\tilde{P}$ is the option impact price and can be either defined by PDE (4.9) or through the reduced to Equation (4.13). Thus,
the cost arising from the strategy $x$ is

$$ \mathcal{C}(x) := \int_{0}^{T} \tilde{P}_{t} \dot{x}_{t} dt $$

The agent's objective is then to minimize a certain objective function, which takes into account his risk aversion, and may involve both cost and risk terms, over the class of admissible trading strategies $x$ with side conditions $x_0 = X$ and $x_T = 0$. This is known as the *optimal trade execution problem*.

In this paper we will treat two cases:

- The expected cost $\mathbb{E}[\mathcal{C}(x)]$

- The mean-variance case $\mathbb{E}[\mathcal{C}(x)] + \lambda \mathbb{V}ar[\mathcal{C}(x)]$

Except from Section 4.6, we will develop the framework under the Black & Scholes case. $\tilde{P}_t$ will be defined by Equation (4.10), and permanent impact will be excluded, i.e $\tilde{\gamma} = 0$.

### 4.3.2 The Black-Scholes framework under a temporary market impact

The effective price $\tilde{P}_t$ is given by Equation (4.13) with $\tilde{\gamma} = 0$ is reduced to:

$$ \tilde{P}_t = P_t + \frac{1}{2} \tilde{\eta} \dot{x}_t \sigma S_t^2 (\hat{T} - t)^{3/2} \Gamma(t, S_t). \quad (4.14) $$

We can rewrite the cost function as the following:

$$ \mathcal{C}(x) = \int_{0}^{T} P_{t} \dot{x}_{t} dt + \frac{1}{2} \tilde{\eta} \int_{0}^{T} \dot{x}_{t}^{2} \sigma S_{t}^{2} (\hat{T} - t)^{3/2} \Gamma(t, S_{t}) dt. $$

Using a simple integration by part and Ito's formula, the cost arising from the strategy $x$ becomes:

$$ \mathcal{C}(x) = -XP_0 - \int_0^T \sigma x_t S_t \Delta(t, S_t) dW_t + \frac{1}{2} \tilde{\eta} \sigma \int_0^T \dot{x}_t^2 S_t^2 (\hat{T} - t)^{3/2} \Gamma(t, S_t) dt, \quad (4.15) $$

where $\Delta$ is the Black-Scholes delta of the option calculated on $\sigma$.

The expected cost of strategy $x$ is then:

$$ \mathbb{E}[\mathcal{C}(x)] = -XP_0 + \frac{1}{2} \tilde{\eta} \mathbb{E} \left[ \int_0^T \dot{x}_t^2 S_t^2 (\hat{T} - t)^{3/2} \Gamma(t, S_t) dt \right]. \quad (4.16) $$

**Theorem 1.** The optimal strategy $x^*$ resulting in minimizing the expected cost under the Black & Scholes framework given by Equation (4.16) is characterized by:

$$ \begin{aligned} \dot{x}^*(t) &= \frac{K_1}{(\hat{T} - t)^{3/2}} \\ x^*(t) &= \frac{K_1}{(\hat{T} - t)^{1/2}} + K_2 \end{aligned} $$

$$ \text{where } K_1 = \frac{X}{2(\hat{T}^{-1/2} - (\hat{T}-T)^{-1/2})} \quad \text{and} \quad K_2 = -2K_1(\hat{T}-T)^{-1/2}. $$
**Démonstration.** Under the Black-Scholes framework, $S_t^2 \Gamma(t, S_t)$ is a martingale which allows to write

$$ \min_x \mathbb{E}[\mathcal{C}(x)] = \frac{1}{2} \tilde{\eta} S_0^2 \Gamma(0, S_0) \min_x \mathbb{E} \left[ \int_0^T (\hat{T} - t)^{3/2} \dot{x}_t^2 dt \right]. $$

Assuming the strategy $x$ is deterministic, the problem is reduced to solving

$$ \min_x \int_0^T (\hat{T} - t)^{3/2} \dot{x}(t)^2 dt. \qquad (4.17) $$

Using the calculus of variations, we find that the solution to (4.17) verifies the Euler-Lagrange equation

$$ \frac{d}{dt} (2(\hat{T} - t)^{3/2} \dot{x}(t)) = 0, $$

along with the boundary conditions

$$ x(0) = X, \quad x(T) = 0. $$

Which gives

$$ \begin{aligned} \dot{x}^*(t) &= \frac{K_1}{2(\hat{T} - t)^{3/2}} \\ x^*(t) &= \frac{K_1}{(\hat{T} - t)^{1/2}} + K_2, \end{aligned} $$

where $K_1 = \frac{X}{2(\hat{T}^{-1/2} - (\hat{T}-T)^{-1/2})}$ and $K_2 = -2K_1(\hat{T}-T)^{-1/2}$.

Moreover, using the Cauchy-Schwarz inequality

$$ \int_0^T f^2(t) dt \int_0^T g^2(t) dt \geq \left( \int_0^T f(t) g(t) dt \right)^2, $$

where $f(t) = \frac{1}{(\hat{T}-t)^{3/4}}$ and $g(t) = (\hat{T}-t)^{3/4} \dot{x}_t$, we have

$$ \mathbb{E} \left[ \int_0^T (\hat{T} - t)^{3/2} \dot{x}_t^2 dt \right] \geq \frac{X^2}{\int_0^T (\hat{T} - t)^{-3/2} dt}. $$

And the equality holds for $(\dot{x}^*, x^*)$.

□

**Remarks 3. :**

- We recall that the expected cost optimal strategy for the equity case is characterized by having a constant trading rate $\dot{x}_t^* = -\frac{X}{\hat{T}}$, as shown in [27] in a discrete-time setting. In the option framework under the impact function we select, the trading speed is an increasing convex function of time.

- The expected cost under the Black & Scholes framework is the only case where a closed solution can be found. This strategy is illustrated in Figure 4.3 for $\hat{T} = 1$, $T = 0.5$ and $X = -1$. We see that both the trading rate and the inventory are increasing and convex. It is to mention that the optimal strategy arising from the expected cost framework does not depend on the underlying asset price.
<figure><img src="image_49.png" /><figcaption>FIGURE 4.3 – The trading strategy for an expected cost framework under the Black-Scholes setting.</figcaption></figure>

In the following sections, we develop the optimal execution framework under a risk/reward criterion under the Black & Scholes framework and the impact price defined by Equation (4.10).

## 4.4 Adding the agent risk aversion : A mean-variance framework

The expected cost is usually used for an agent who does not monitor the risk of his strategy. Investors, however, usually takes into account their risk aversion, through utility functions or using risk/reward criterion such as the mean-variance. The literature on these problems is rich for optimal execution of a book of equity shares. For example Almgren and Chriss [7], and Forsyth [75], study the mean-variance optimal execution problem. In [90], Gatheral and Schied take the time-average value-at-risk associated with the P&L of the position, while Forsyth and co-authors, see [76], use a quadratic variance as a risk criterion. In this paper, we focus on the mean-variance criterion in light of [6, 7]. The mean-variance of the cost of trading is defined by :

$$ \mathbb{E}[\mathcal{C}(x)] + \tilde{\lambda}\mathbb{V}[\mathcal{C}(x)], $$

where $\tilde{\lambda} > 0$ is the variance penalty. The choice of this coefficient can not be explained in terms of fundamental investment preferences. The value is chosen in order to obtain solutions that bring out a certain meaning to the optimization problem.

The variance of the cost function term can be written as the following

$$
\begin{align*}
\mathbb{V}[\mathcal{C}(x)] &= \mathbb{E}\left[\left(\int_0^T \tilde{P}_t \dot{x}_t dt - \mathbb{E}\left[\int_0^T \tilde{P}_t \dot{x}_t dt\right]\right)^2\right] \\
&= \mathbb{E}\left[\left(\int_0^T \frac{1}{2} x_t (\tilde{\sigma}_t^2 - \sigma^2) S_t^2 \partial_{SS} \tilde{P}(t, S_t) dt - \int_0^T x_t \sigma S_t \partial_S \tilde{P}(t, S_t) dW_t\right)^2\right] \\
&= \mathbb{E}\left[\int_0^T x_t^2 \sigma^2 S_t^2 \partial_S \tilde{P}^2(t, S_t) dt\right] \\
&\quad + \text{\{terms arising from uncertainty in the drift part\}}.
\end{align*}
$$

The exact expression of the variance of $\mathcal{C}(x)$ is complicated since all terms are random. A reasonable assumption is that the largest source of uncertainty arises from the stochastic integral part.
In [121, 153], the authors explore this approximation. They argued that the terms in the drift part are small compared with market dynamics. In fact, we verify numerically that such approximation makes sense in our case. We consider a time discretization of the trading interval $[0, T]$ and compare the "true" variance to the approximation by taking the stochastic integral part. The error is indeed small due to the coefficient $\bar{\eta}^2$ in the drift part.

In this section we are interested in the price impact formulation of Equation (4.10) with temporary impact only. That is, we can easily deduce that the mean-variance objective function can be approximated as the following:

$$ \mathbb{E}[\mathcal{C}(x)] + \lambda \mathbb{E}[\mathcal{C}(x)] \approx \mathbb{E}\left[\int_{0}^{T} \frac{1}{2}\bar{\eta}\sigma\dot{x}_{t}^{2}S_{t}^{2}(\hat{T} - t)^{3/2}\Gamma(t, S_{t})dt + \bar{\lambda}\int_{0}^{T} x_{t}^{2}\sigma^{2}S_{t}^{2}\Delta^{2}(t, S_{t})dt\right], $$

where $\Delta(t, S_t)$ is the Black-Scholes delta.

We define $\mathcal{X}$ (T,X) the set of all adapted and absolutely continuous strategies that satisfy the boundary conditions $x_0 = X < 0$, $x_T = 0$, the integrability conditions $\mathbb{E}[\int_0^T \dot{x}_t^2 S_t^2 (\hat{T} - t)^{3/2} \Gamma(t, S_t) dt] < \infty$ and $\mathbb{E}[\int_0^T x_t^2 S_t^2 \Delta^2(t, S_t) dt] < \infty$, and consider the following minimization problem

$$ U(0, S_0, X) := \tag{4.18} $$

$$ \inf_{x \in \mathcal{X}(T,X)} \mathbb{E} \left[ \int_0^T \left\{ \dot{x}_t^2 S_t^2 (\hat{T} - t)^{3/2} \Gamma(t, S_t) + \lambda x_t^2 \sigma^2 (S_t) S_t^2 \Delta^2(t, S_t) \right\} dt \right], \tag{4.19} $$

where $\lambda := 2\bar{\lambda}/\sigma\bar{\eta}$.

This minimization problem does not admit a closed-form solution. In order to solve it, a proper stochastic dynamic programming framework needs to be set up.

**Proposition 3.** We parameterize strategies $x$ by their speed of trading and define the control $\alpha$ such that $\alpha_t := -\dot{x}_t$. We introduce $\mathcal{A}$ (T,X) the class of all progressively measurable processes $(\alpha_t)_{0 \le t \le T}$, for which the parameterized strategy $x^\alpha$ defined by

$$ x_t^\alpha := X - \int_0^t \alpha_s ds, \quad 0 \le t \le T, $$

belongs to the set $\mathcal{X}$ (T,X). We restrict our framework to Markovian controls and thus, solving such optimal stochastic control problem at time 0 is brought to a more general case where the agent starts buying at any arbitrary time $t \in [0, T]$ with an initial quantity $x$ without losing the optimality.

$$ \alpha_t = \alpha(t, S_t, x_t) $$

**Remark 8.** The agent's optimal trading speed $\alpha_t$ at time $t$ is completely determined by the current state $(t, S_t, x_t)$, i.e time $t$, current stock price $S_t$ and current quantity $x_t$.

For a given strategy $\alpha(.,,.)$, the value function $U(t, S, x)$ is defined as

$$ U(t, S, x) = \inf_{\alpha \in \mathcal{A}(T,X)} \mathbb{E}_t \left[ \int_t^T \left\{ \alpha_u^2 S_u^2 (\hat{T} - u)^{3/2} \Gamma(u, S_u) + \lambda \sigma^2 (x_u^\alpha)^2 S_u^2 \Delta^2(u, S_u) \right\} du \right], $$

where $\mathbb{E}_t$ is the expectation conditional to $S_t = s$ and $x_t^\alpha = x$. Note that this problem fits into the finite-fuel framework because of the state dependence of the class $\mathcal{A}$ (T,X) of admissible control processes.
**Remarks 4.** (i) We can rewrite the state variables $S_t$ and $x_t^\alpha$ in terms of stochastic differential equations (SDE):

$$dS_t = \sigma S_t dW_t$$
$$dx_t^\alpha = -\alpha_t dt$$

where the initial conditions are $S_0 = s_0$ and $x_0 = X$.

(ii) Using the standard procedure of deriving the Hamilton-Jacobi-Bellman (HJB) equation in stochastic control problems, see [157], $U(t, S, x)$ should satisfy the following HJB equation

$$\partial_t U + \frac{1}{2}\sigma^2 S^2 \partial_{SS} U + \bar{\lambda}x^2 S^2 \Delta^2(t, S) + \inf_{\alpha \in \mathbb{R}} \{\alpha^2 S^2 (\hat{T} - t)^{3/2} \Gamma(t, S) - \alpha \partial_x U\} = 0. \quad (4.20)$$

(iii) Assuming $U$ is sufficiently smooth, Bellman's principle of optimality for stochastic minimization problems suggests that the value function of the dynamic programming is always a submartingale and a martingale under the optimal strategy.

(iv) The so-called finite-fuel constraint required from strategies in $\mathcal{A}(T,X)$ (i.e $\int_0^T \alpha_t dt = X$), suggests the value function $U$ should satisfy a singular terminal condition of the form

$$\lim_{t \to T} U(t, S, x) = \begin{cases} 0 & \text{if } x=0 \\ +\infty & \text{if } x \neq 0. \end{cases} \quad (4.21)$$

The intuition of Equation (4.21) is that a state with a non zero option position with no time left for its liquidation means that the liquidation task has not been performed, and so it should receive an infinite penalty.

In what follows, we will substitute the infinite penalty problem with a finite terminal condition. We index the value function by $\epsilon$ such that

$$U_{\epsilon}(t, s, x) = \inf_{\alpha \in \mathcal{A}(T,X)} \mathbb{E}_{t} \left[ \int_{t}^{T} \left\{ \alpha_{u}^{2} S_{u}^{2} (\hat{T} - u)^{3/2} \Gamma(u, S_{u}) + \lambda \sigma^{2} (x_{u}^{\alpha})^{2} S_{u}^{2} \Delta^{2}(u, S_{u}) \right\} du + \frac{1}{\epsilon} \psi(x_{T}^{\alpha}) \right]. \quad (4.22)$$

With terminal condition

$$U_{\epsilon}(T, s, x) = \frac{1}{\epsilon} \psi(x) \begin{cases} 0 & \text{if } x=0 \\ \gg 1 & \text{if } x \neq 0. \end{cases} \quad (4.23)$$

**Remarks 5.** (i) The function $\psi$ should equal to 0 when $x=0$ and at the same time we need to be able to penalize a final state with remaining inventory through $\epsilon$. A proper $\psi$ will be more convenient than another one. In our case, the dependence on the state variable $x$ suggests $\psi(x) = x^2$. This ansatz is used in [104] for example.

(ii) The solution to (4.22)-(4.23) is $\epsilon$-dependent. For different values $\epsilon$ and $\epsilon'$ we find different solutions $U_\epsilon$ and $U_{\epsilon'}$. However, we can construct a time series $(U_{\epsilon_n})_n$ that converges to the original problem when $\epsilon_n \xrightarrow{n \to +\infty} +\infty$.

(iii) Finally, we can easily verify that $U_\epsilon$ satisfies the HJB equation (4.20) with terminal condition (4.23).

In the following, we seek to solve the dynamic programming problem (4.22) with terminal condition (4.23). Our first main result is the following:

**Theorem 2.** Let $U_\epsilon^*$ be a regular function which solves the PDE:

$$\left\{ \begin{array}{l} \partial_t U_\epsilon^* + \frac{1}{2}\sigma^2 S^2 \partial_{SS} U_\epsilon^* + \lambda x^2 \sigma^2 S^2 \Delta^2(t, S) - \frac{(\partial_x U_\epsilon^*)^2}{4(\hat{T}-t)^{3/2} I(t,S)} = 0 \\ U_\epsilon^*(T, S_T, x_T) = \frac{1}{\epsilon} \psi(x_T^\alpha). \end{array} \right. \quad (4.24)$$
Then $U_{\varepsilon}^*$ is the unique solution to the optimal execution problem (4.22). Moreover, the optimal execution rate $\alpha_t^* = -\dot{x}_t^*$ is such that:

$$ \alpha_t^* = \frac{\partial_x U_{\varepsilon}^*(t, S_t, x_t^*)}{4(\hat{T} - t)^{3/2} S_t^2 \Gamma(t, S_t)}. \quad (4.25) $$

The proof of this theorem uses the following lemma:

**Lemma 1.** Let $U_{\varepsilon}^*$ be a regular function that solves PDE (4.24). Then $U_{\varepsilon}^*$ is a classical solution of the minimization problem (4.20), (4.23), and the optimal execution speed $\alpha_t^* = -\dot{x}_t^*$ is defined by:

$$ \begin{aligned} \alpha_t^* &= \frac{\partial_x U^*(t, S_t, x_t^*)}{4(\hat{T} - t)^{3/2} S_t^2 \Gamma(t, S_t)} \quad (4.26) \\ &= \operatorname{argmin} \left\{ \alpha^2 (\hat{T} - t)^{3/2} S_t^2 \Gamma(t, S_t) - \alpha \partial_x U^*(t, S_t, x_t^*) \right\}. \end{aligned} $$

*Démonstration.* The proof to the lemma is straightforward and uses the fact that the function $h(\alpha) = \alpha^2(\hat{T}-t)^{3/2}S_t^2\Gamma(t,S_t)-\alpha\partial_x U_{\varepsilon}^*(t,S_t,x_t^*)$ is convex and attains its minimum for $h'(\alpha) = 0$. This gives (4.26) and

$$ h(\alpha^*) = - \frac{\left(\partial_x U^*(t, S_t, x_t^*)\right)^2}{4(\hat{T} - t)^{3/2} S_t^2 \Gamma(t, S_t)}. $$

Thus, we obtain the lemma through substituting $\alpha^*$ back into Equation (4.20). □

**Proof of Theorem 2:** Let $\alpha \in \mathcal{A}(T,X)$ be an arbitrary control process, and define the stopping time

$$ \theta_n := T \wedge \inf \left\{ s > t : \int_t^s \alpha_u du \ge n \right\} \quad (4.27) $$

We first need to show that

$$ U_{\varepsilon}^{*}(t, S, x) \le \mathbb{E} \left[ \int_{t}^{T} \left\{ \alpha_{u}^{2}(\hat{T} - u)^{3/2} S^{2} \Gamma(u, S_{u}) + \lambda \sigma^{2} (x_{u}^{\alpha})^{2} S_{u}^{2} \Delta^{2}(u, S_{u}) \right\} du + \frac{1}{\varepsilon} \psi(x_{T}^{\alpha}) \right], $$

where $U_{\varepsilon}^*$ is the solution of the PDE (4.24).

First, we have that the right hand side of the inequality is finite. This is justified by the admissibility of $\alpha$ which includes the integrability conditions $\mathbb{E}[\int_t^T \alpha_u^2(\hat{T}-u)^{3/2}S_u^2\Gamma(u,S_u)du] < \infty$ and $\mathbb{E}[\int_t^T (x_u^\alpha)^2 S_u^2 \Delta^2(u,S_u)du] < \infty$. In addition, the fuel constraint $\int_0^T \alpha_s ds = X$, and Cauchy-Schwarz inequality imply

$$ |x_t^\alpha| = \left| \int_t^T \alpha_s ds \right| \le \sqrt{(T-t) \int_0^T \alpha_s^2 ds} \in L^2(\mathbb{P}). $$

By Ito's formula, we have

$$ \begin{aligned} U_{\varepsilon}^{*}(t, S, x) = U_{\varepsilon}^{*}(\theta_n, S_{\theta_n}, x_{\theta_n}) & - \int_{t}^{\theta_n} \left( \partial_t + \frac{1}{2} \sigma^2 S_t^2 \partial_{SS} - \alpha_r \partial_{X} \right) U_{\varepsilon}^{*}(r, S_r, x_r) dr \\ & - \int_{t}^{\theta_n} \sigma S_r \partial_{S} U_{\varepsilon}^{*}(r, S_r, x_r) dW_r. \end{aligned} $$

Observe that

$$ \begin{aligned} & \partial_t U_{\varepsilon}^* + \frac{1}{2} \sigma^2 S_r^2 \partial_{SS} U_{\varepsilon}^* - \alpha_r \partial_{X} U_{\varepsilon}^* + \alpha_r^2 (\hat{T} - r) S^2 \Gamma(r, S_r) + \lambda \sigma^2 (x_r^\alpha)^2 S_r^2 \Delta^2(r, S_r) \\ & \ge \partial_t U_{\varepsilon}^* + \frac{1}{2} \sigma^2 S^2 \partial_{SS} U_{\varepsilon}^* + \lambda \sigma^2 (x_r^\alpha)^2 S_r^2 \Delta^2(r, S_r) + \inf_{\alpha \in \mathcal{A}(X,T)} \{\alpha^2 (\hat{T} - r) S^2 \Gamma(r, S_r) - \alpha \partial_{X} U_{\varepsilon}^*\}, \end{aligned} $$
and that the integrand in the stochastic integral is bounded on $[t, \theta_n]$, a consequence of the regularity of $U_{\varepsilon}^*$. It follows that:

$$U_{\varepsilon}^{*}(t, S, x) \leq \mathbb{E}[U_{\varepsilon}^{*}(\theta_{n}, S_{\theta_{n}}, x_{\theta_{n}})] + \mathbb{E}\left[\int_{t}^{\theta_{n}} \{\alpha_{r}^{2}(\hat{T} - r)S^{2}\Gamma(r, S_{r}) + \lambda\sigma^{2}(x_{r}^{\alpha})^{2}\Delta^{2}(r, S_{r})\}dr\right].$$

We now take the limit as $n$ increases to infinity. Since $\theta_n \to T$ a.s we have

$$\lim_{n \to \infty} \mathbb{E} \left[ \int_{t}^{\theta_n} \{ \alpha_r^2 (\hat{T} - r) S^2 \Gamma(r, S_r) + \lambda \sigma^2 (x_r^\alpha)^2 \Delta^2 (r, S_r) \} dr \right] = \mathbb{E} \left[ \int_{t}^{T} \{ \alpha_r^2 (\hat{T} - r) S^2 \Gamma(r, S_r) + \lambda \sigma^2 (x_r^\alpha)^2 \Delta^2 (r, S_r) \} dr \right].$$

On the other hand

$$\begin{aligned} |U_{\varepsilon}^{*}(\theta_{n}, S_{\theta_{n}}, x_{\theta_{n}})| & \leq \frac{1}{\varepsilon}\psi(x_{T}^{\alpha}) + K_{1}\max_{t \leq r \leq T}(S_{r}^{2}) + K_{2}\max_{t \leq r \leq T}(S_{r}^{2})\int_{t}^{T}(x_{r}^{\alpha})^{2}dr + K_{3} \\ & \leq \frac{1}{\varepsilon}\psi(x_{T}^{\alpha}) + K_{1}\max_{t \leq r \leq T}(S_{r}^{2}) + \frac{1}{2}K_{2}\max_{t \leq r \leq T}(S_{r}^{2})T^{2}\int_{0}^{T}\alpha^{2}dr + K_{3}. \end{aligned}$$

It follows from the dominated convergence

$$U_{\varepsilon}^{*}(t, S, x) \leq \mathbb{E} \left[ \int_{t}^{T} \left\{ \alpha_{u}^{2}(\hat{T} - u)S^{2}\Gamma(u, S_{u}) + \lambda\sigma_{0}^{2}(x_{u}^{\alpha})^{2}S_{u}^{2}\Delta^{2}(u, S_{u}) \right\} du + \frac{1}{\varepsilon}\psi(x_{T}^{\alpha}) \right].$$

The control given by 4.25 is well defined ($T<\hat{T}$) and the solution $U_{\varepsilon}$ regular.

**Remark 9.** Setting the control variable as the speed of trading $\dot{x}$ is very common. In [104] for example, the authors use this parametrization to obtain the HJB equation. The problem is then reduced to one dimension using an Ansatz. An alternate solution is to adopt an exponential growth parametrization. This way, reducing the dimension is straightforward.

### 4.4.1 Deriving the HJB with multiplicative state variable in the Black-Scholes case

Let us consider the controlled state variable $x_t^K$ and the control $\kappa$ such that

$$dx_t^K = -\kappa_t x_t^K dt$$

where $\kappa_t > 0$, $x_0 = X$ and $x_t$ increasing and bounded by 0. And let us define $\mathcal{K}(T, X)$ the set of admissible control processes $\kappa$ such that $x$ belongs to $\mathcal{X}(T, X)$.

**Remark 10.** Ideally, $x_t$ vanishes as $t$ goes to $T$ to verify the finite fuel constraint of total acquisition (i.e liquidation of a short position) at $T$. This is only possible if $\kappa$ is infinite at a certain time. As in the previous case, the limit $\varepsilon \to 0$ forbids to trade a large quantity at the end time by imposing a penalty for the final state. Substituting the infinite limit by a large penalty on the value function as $t \to T$ allows to have a regular control variable which can be solved numerically. The remaining inventory will be acquired as a last additional trade at $t=T$.

Using the multiplicative parameterization we have:

$$\begin{aligned} & \mathbb{E}_t \left[ \int_t^T \left\{ \kappa_u^2 (x_u^\kappa)^2 (\hat{T} - u)^{3/2} S_u^2 \Gamma(u, S_u) + \lambda (x_u^\kappa)^2 \sigma^2 S_u^2 \Delta^2 (u, S_u) \right\} du \right] \\ & = x^2 \mathbb{E}_t \left[ \int_t^T e^{-\int_t^u 2\kappa_s ds} \left\{ \kappa_u^2 (\hat{T} - u)^{3/2} S_u^2 \Gamma(u, S_u) + \lambda \sigma^2 S_u^2 \Delta^2 (u, S_u) \right\} du \right]. \end{aligned} \quad (4.28)$$
Moreover, by taking the terminal condition $U_{\varepsilon}(T, S_T, x_T) = \frac{1}{\varepsilon}(x_T^{\kappa})^2 = x^2 \frac{1}{\varepsilon} e^{-\int_t^T 2\kappa_s ds}$ we have:

$$U_{\varepsilon}(t, s, x) := x^2 u_{\varepsilon}(t, s)$$

where $u_{\varepsilon}$ is a reduced value function.

**Definition 1.** The reduced value function $u_{\varepsilon}$ is defined by:

$$u_{\varepsilon}(t, S) = \inf_{\kappa \in \mathcal{K}} \mathbb{E}_t \left[ \int_t^T e^{-\int_t^u 2\kappa_s ds} \{\kappa_u^2 (\hat{T} - u)^{3/2} S_u^2 \Gamma(u, S_u) + \lambda \sigma^2 S_u^2 \Delta^2(u, S_u)\} du + \frac{1}{\varepsilon} e^{-\int_t^T \kappa_s ds} \right]. \quad (4.29)$$

**Proposition 4.** By means of Ito's formula, $u_{\varepsilon}$ verifies the HJB equation:

$$\partial_t u_{\varepsilon} + \frac{1}{2} \sigma^2 S^2 \partial_{SS} u_{\varepsilon} + \inf_{\kappa} \{\kappa^2 (\hat{T} - t)^{3/2} S^2 \Gamma(t, S) - 2\kappa u_{\varepsilon}\} + \lambda \sigma^2 S^2 \Delta^2(t, S) = 0 \quad \text{and} \quad u_{\varepsilon}(T, s) = \frac{1}{\varepsilon}. \quad (4.30)$$

Besides, $h(\kappa) = \kappa^2 (\hat{T} - t)^{3/2} S^2 \Gamma(t, S) - 2\kappa u_{\varepsilon}(t, S)$ attains its minimum for $h'(\kappa) = 0$ and

$$\kappa^*(t, S) = \frac{u_{\varepsilon}(t, S)}{(\hat{T} - u)^{3/2} S^2 \Gamma(t, S)} \quad \text{and} \quad h(\kappa^*) = - \frac{u_{\varepsilon}^2}{(\hat{T} - u)^{3/2} S^2 \Gamma(t, S)}.$$

By injecting the previous expression into the HJB equation (4.30), we deduce the PDE for $u_{\varepsilon}$.

**Theorem 3.** Let $u_{\varepsilon}^*$ be a regular function verifying the following PDE

$$\left\{ \begin{array}{l} \partial_t u_{\varepsilon}^* + \frac{1}{2} \sigma^2 S^2 \partial_{SS} u_{\varepsilon}^* + \lambda \sigma^2 S^2 \Delta^2(t, S) - \frac{1}{(\hat{T} - t)^{3/2} S^2 \Gamma(t, S)} u_{\varepsilon}^2 = 0 \\ u_{\varepsilon}^*(T, s) = \frac{1}{\varepsilon}. \end{array} \right. \quad (4.31)$$

Then $u_{\varepsilon}^*$ is the unique solution to the reduced optimization problem (4.29). The optimal trading rate $\kappa_t^*$ is defined by:

$$\kappa^*(t, S) = \frac{u_{\varepsilon}(t, S)}{(\hat{T} - t)^{3/2} S^2 \Gamma(t, S)}.$$

**Remarks 6.** (i) The proof of Theorem 3 uses the same verification arguments as developed for Theorem 2. (ii) The multiplicative form gives an interesting interpretation for the control variable. By linearizing the quadratic term in PDE (4.31), we can see $\frac{1}{(\hat{T} - t)^{3/2} S^2 \Gamma(t, S)}$ $u_{\varepsilon}$ as a discount factor of the value function. In the numerical scheme, this term will be taken explicitly in order to solve the PDE. (iii) There is a feedback between the value function $u_{\varepsilon}$ and the control $\kappa$. (iv) Using the analogy to equity markets, the denominator $(\hat{T} - t)^{3/2} S^2 \Gamma(t, S)$ in the optimal trading rate can be interpreted as a market depth. In our case, the depth is a function of time and asset price.

### 4.4.2 Localisation and boundary conditions

The original reduced problem (4.29) is posed on the domain $(t, s) \in [0, T] \times [0, \infty]$. For computational purposes, and because asset prices are finite, one needs to localize this domain to $[0, T] \times [0, S_{\text{max}}]$. Thus, we need to add the following complementary conditions:

— When $S = 0$, the put price tends to the strike : $P(t, S) \approx K$, $S^2 \Gamma(t, S) \approx 0$ and $S^2 \Delta^2(t, S) \approx 0$. We simply need to solve

$$\partial_t u_{\varepsilon} + \inf_{\kappa}\{-2\kappa u_{\varepsilon}\} = 0 \quad \text{and} \quad u_{\varepsilon}(T) = \frac{1}{\varepsilon}.$$

This limit condition is a singular control problem which can also be expressed as a variational equation

$$\min \{\partial_t u_{\varepsilon}, u_{\varepsilon}\} = 0 \quad \text{and} \quad u_{\varepsilon}(T) = \frac{1}{\varepsilon}.$$
— When $S = S_{\max} \gg K$, the put price becomes insignificant. We have $S^2\Gamma(t, S) \approx 0$ and $S^2\Delta^2(t, S) \approx 0$, which makes both the impact and the risk aversion term vanish. The optimization reduces to solving

$$ \partial_t u_{\epsilon} + \frac{1}{2}\sigma^2 S^2 \partial_{SS} u_{\epsilon} + \inf_{\kappa}\{-2\kappa u_{\epsilon}\} = 0 \text{ and } u_{\epsilon}(T) = \frac{1}{\epsilon}. $$

The variational inequality arising from this condition is:

$$ \min \left\{ \partial_t u_{\epsilon} + \frac{1}{2} \sigma^2 S^2 \partial_{SS} u_{\epsilon}, u_{\epsilon} \right\} = 0 \text{ and } u_{\epsilon}(T) = \frac{1}{\epsilon}. $$

By taking these approximations on the boundaries, we are able to solve the problem numerically using finite differences methods as we will see in the next section.

# 4.5 Numerical solution and results

## 4.5.1 A finite differences scheme

The PDE (4.31) with additional boundary conditions of Section 4.4.2 is a Riccati equation. In [146], the authors solve the Algmren-Chriss optimal execution problem in terms of the Riccati equation using two methods. They find the exact solution using a time "reparametrization". Unfortunately, our case does not allow for a closed-form solution. Thus, finite differences methods is a convenient way to solve the problem.

Restricting the problem to the domain $[0, T] \times [0, S_{\max}]$, we can use the argument of Lipschitz functions on a compact to solve the semilinear PDE (4.31) numerically. We discretize time and space and define $\tau = T - t$ as the time to the strategy end time. We fix $\epsilon$ and denote by $u_j^n$ the numerical approximation to $u_{\epsilon}(n\Delta\tau, j\Delta S)$, where $j = 1, 2, \dots, J-1$ is the space grid index, and $n = 1, 2, \dots, N-1$ the time grid index taken backwards. $\Delta S$ is the spacial step size and $\Delta\tau = \Delta t$ the time step.

Let $\mathcal{L}_h u_j^n$ denote the spatial discretization of the differential term $\mathcal{L} u = \frac{1}{2}\sigma^2 \partial_{SS} u$, where:

$$ \mathcal{L}_h u_j^n = \frac{1}{2} \sigma^2 S_j^2 \frac{u_{j+1}^n - u_j^n + u_{j-1}^n}{\Delta S^2} \quad (4.32) $$

We linearize the quadratic term by decomposing it into the product of an explicit and an implicit form. Thus, the general family of the two-level implicit schemes for solving the equation is given by:

$$ \frac{u_j^{n+1} - u_j^n}{\Delta\tau} = \theta \mathcal{L}_h u_j^{n+1} + (1-\theta) \mathcal{L}_h u_j^n + \frac{1}{(\hat{T} - t_n)^{3/2} S_j^2 \Gamma_j^n} u_j^n u_j^{n+1} + \lambda \sigma^2 (S_j)^2 (\Delta_j^n)^2, $$

where $\Gamma_j^n = \Gamma(n\Delta\tau, j\Delta S)$ and $\Delta_j^n = \Delta(n\Delta\tau, j\Delta S)$.

For terminal condition $n=0$ and boundary conditions $j=0$ and $j=J$ we have the following:

— $u_j^0 = \frac{1}{\epsilon}$ where $\frac{1}{\epsilon} \gg 1$ which translates the penalty related to the finite fuel constraint.

— $u_0^n$ verifies:

$$ \min \left( - \frac{u_0^{n+1} - u_0^n}{\Delta\tau}, u_{\epsilon 0}^n \right) = 0. $$
<table><caption>TABLEAU 4.1 – Parameters for buying options under market impact over 1 month horizon</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>σ</td><td>30%</td></tr><tr><td>T (the strategy horizon)</td><td>1/12 (years)</td></tr><tr><td>Τ̂ (the option maturity)</td><td>1 (years)</td></tr><tr><td>μ</td><td>0</td></tr><tr><td>r</td><td>0</td></tr><tr><td>S<sub>0</sub></td><td>1</td></tr><tr><td>K</td><td>S<sub>0</sub></td></tr><tr><td>Action</td><td>Buy</td></tr><tr><td>x<sub>0</sub></td><td>-1</td></tr><tr><td>η̄</td><td>0.05</td></tr><tr><td>Trading frequency</td><td>4 trades per day</td></tr><tr><td>λ</td><td>0, 1, 10, 100</td></tr></tbody></table>

— And $u_j^n$ verifies :

$$ \min \left( -\frac{u_j^{n+1} - u_j^n}{\Delta \tau} - \theta \mathcal{L}_h u_j^{n+1} - (1-\theta) \mathcal{L}_h u_j^n, u_j^n \right) = 0. $$

The inventory $x_n$ at $\tau_n$ is expressed by

$$ x_n^j = x_{n+1}^j e^{\kappa_n^j \Delta \tau} $$

where $x_N = X$, $x_0 = 0$ and the control is given by

$$ \kappa_{t_n}^j = \frac{u_n^j}{(\hat{T} - t_n)^{3/2} S_j^2 \Gamma_j^n} $$

Thus, the quantity traded at $\tau_n$ at the asset price node $S_j$ is

$$ x_n^j - x_{n+1}^j = x_{n+1}^j (e^{\kappa_n^j \Delta \tau} - 1). $$

### 4.5.2 Results

Usually end-users have a net long positions in OTM or ATM puts. These positions are explained by the fact that end-users suffer from "crashophobia" as explained in [147]. In our numerical experiment, we present results for a long position on ATM put options. The moneyness is fixed w.r.t the asset price at time 0. The strike K is set to be equal $S_0$ and remains the same until the strategy ending date. The parameters for this example are given in Table 4.1.

Figure 4.4 gives the optimal execution strategy through the rate of trading $\kappa$ as a function of the underlying price $S$ and time $t$. The strategy does not depend on the trader position. However, as time increases the trading rate increases as well. The dependence of $\kappa$ on the underlying price direction is barely notable. First, the case $\lambda = 0$ corresponds to the expected cost which was found to give a strategy that does not depend on the underlying price. The trading speed $\dot{x}_t$ of the strategy is increasing and convex in time, with contrast to the equity case where it was found to be constant. The mean-variance adds some dependency on the asset level but the surface representation does not allow to see it.
<figure><img src="image_50.png" /><figcaption>FIGURE 4.4 – The rate of trading $κ$ as a function of the underlying price $S$ and time $t$ for different values of $λ$ ($λ = 0$ top left, $λ = 1$ top right, $λ = 10$ bottom left, $λ = 100$ bottom right). The strike $K = S_0$ is fixed at time 0.</figcaption></figure>

Figure 4.4 is misleading in a way; when the agent runs the strategy, the inventory and fundamental price will all evolve, hence this representation is to consider carefully. To gain additional insight into the dynamic behavior, we plot in Figures 4.5, 4.6 and 4.7 4 paths of the underlying price together with the rate of trading $κ$, the inventory $x$ and quantity to be traded $\Delta x$. We can see clearly that adding the variance pushes the agent to adapt the strategy to the underlying level (Figure 4.6). Furthermore, when as the risk aversion parameter $λ$ increases, the traded quantity tends to be larger at the beginning (Figure 4.7). Finally, we plot in Figure 4.8 a heat map of the distribution of the trading rate $κ$, trading speed $\dot{x}$, and inventory $x$ over 10,000 simulations. This representation allows to see that the mean-variance with a high risk aversion is most sensitive to price movements. The case $λ = 0$ is the least affected by the spot variation.
<figure><img src="image_51.png" /><figcaption>FIGURE 4.5 – Sample paths of the evolution of the fundamental price, trading rate, inventory and traded quantity throughout the execution for $\lambda = 0$.</figcaption></figure>
<figure><img src="image_52.png" /><figcaption>FIGURE 4.6 – Sample paths of the evolution of the fundamental price, trading rate, inventory and traded quantity throughout the execution for $\lambda = 10$.</figcaption></figure>
<figure><img src="image_53.png" /><figcaption>FIGURE 4.7 – Sample paths of the evolution of the fundamental price, trading rate, inventory and traded quantity throughout the execution for $\lambda = 100$.</figcaption></figure>
<figure><img src="image_54.png" /><figcaption>FIGURE 4.8 – Heat maps showing the density of inventory and trading speed throughout the execution for $\lambda = 0, 10,$ and $100$.</figcaption></figure>
## 4.6 Extension to a local volatility model : A numerical method for the general case

In previous sections we exploited a closed-formula of the impact option price to derive the HJB equation. Under instantaneous impact, we were able to reduce the problem dimension to obtain a one-dimensional PDE which we solved numerically through finite differences methods. This framework was possible by linearizing the impact term in the volatility using a Taylor expansion.

In what follows, we would like to keep the nonlinear formulation of the impact in the asset price through the "enlarged volatility". PDE (4.4) actually provides the price of the option at time $t$ given an inventory $x_t$ and trading speed $\dot{x}_t$ and allows to take a local volatility model as pointed out in Remark 1 and developed in [117]. We will build a numerical scheme in this general case, under a constant elasticity volatility (CEV) model. We rewrite PDE (4.4) of $\tilde{P}$ as the following

$$ \begin{cases} \partial_t \tilde{P}(u, S) + \frac{1}{2} \tilde{\sigma}^2(t, S) S^2 \partial_{SS} \tilde{P}(u, S) = 0, & (u, S) \in [t, \hat{T}[\times]0, \infty] \\ \tilde{P}(\hat{T}, s) = (K - s)^+, & s \in ]0, \infty[, \end{cases} \quad (4.33) $$

where

$$ \tilde{\sigma}^2(t, S) = \sigma^2(S) + (\tilde{\eta}\dot{x}_t + \tilde{\gamma}(x_t - x_0))\sqrt{\hat{T} - t}\sigma(S), $$

and

$$ \sigma(S) = \sigma_0 S^{\beta/2-1}. $$

Recall that we can always obtain the Black-Scholes case by setting $\beta = 2$. In this case, one can directly take the Black-Scholes closed formula with the enlarged volatility.

We are again interested in minimizing the mean-variance of the cost arising from strategy $x$, defined by

$$ \mathcal{C}(x) = \int_{0}^{T} \tilde{P}_{t} \dot{x}_{t} dt. $$

We neglect the drift term in the variance which leads to the following approximation

$$ \forall ar[\mathcal{C}(x)] \approx \mathbb{E} \left[ \int_{0}^{T} x_{t}^{2} \sigma^{2}(S_{t}) S_{t}^{2} \partial_{S} \tilde{P}^{2}(t, S_{t}, x_{t}, \dot{x}_{t}) dt \right]. $$

The mean-variance objective function is thus

$$ \mathbb{E}[\mathcal{C}(x)] + \lambda \forall ar[\mathcal{C}(x)] \approx \mathbb{E} \left[ \int_{0}^{T} \dot{x}_{t} \tilde{P}(t, S_{t}, x_{t}, \dot{x}_{t}) dt + \lambda \int_{0}^{T} x_{t}^{2} \sigma^{2}(S_{t}) S_{t}^{2} \partial_{S} \tilde{P}^{2}(t, S_{t}, x_{t}, \dot{x}_{t}) dt \right]. $$

We would like to find the trading strategy $x$ for the following minimization problem

$$ \inf_{x \in \mathcal{X}(T,X)} \mathbb{E} \left[ \int_{0}^{T} \left\{ \dot{x}_t \tilde{P}(t, S_t, x_t, \dot{x}_t) + \lambda \sigma^2(S_t) x_t^2 S_t^2 \partial_S \tilde{P}^2(t, S_t, x_t, \dot{x}_t) \right\} dt \right], $$

where $\tilde{P}$ follows PDE (4.33). Thus, we develop the dynamic optimization framework for V.
For initial fixed $(t, S, x) = (t, S_t, x_t)$ where $t < T$, trading strategy $\alpha$ such that $dx_t^\alpha = -\alpha_t dt$ and risk aversion $\lambda$, we define the value function $V$ at time $t$

$$ \hat{V}(t, S, x) = \inf_{\alpha} \mathbb{E} \left[ \int_{t}^{T} \left\{ -\alpha_t \tilde{P}(t, S_t, x_t^\alpha, -\alpha_t) + \lambda \sigma^2 (S_t) (x_t^\alpha)^2 S_t^2 \partial_S \tilde{P}^2(t, S_t, x_t^\alpha, -\alpha_t) \right\} dt \right]. \quad (4.34) $$

For $t < T$, let $V = V(\tau = T - t, S, x) = \hat{V}(t, S, x)$. It is easy to find that the optimal control $\alpha^*$ can be obtained by solving the following HJB equation

$$ \partial_{\tau} V = \frac{1}{2} \sigma^2(S) S^2 \partial_{SS} V + \inf_{\alpha} \left\{ -\alpha \partial_x V + -\alpha_t \tilde{P}(t, S, x, -\alpha) + \lambda \sigma^2(S) x^2 S_t^2 \partial_S \tilde{P}^2(t, S, x, -\alpha) \right\}. \quad (4.35) $$

We restrict our variables $\tau$, $S = S(\tau)$ and $x = x(\tau)$ to the domain

$$ \Omega = [0, T] \times [0, S_{\max}] \times [X, 0]. $$

And set the initial condition and boundary conditions for $V$ as suggested in Section 4.4. That is:

— For $\tau = 0$

$$ V(\tau = 0, S, x) \begin{cases} 0 & \text{if } x = 0 \\ \gg 1 & \text{if } x \neq 0. \end{cases} \quad (4.36) $$

— For $S = 0$ we simply solve

$$ \partial_{\tau} V = \inf_{\alpha} \{-\alpha \partial_x V\}. $$

— For $S = S_{\max}$ the function $g$ vanishes as the put price and delta tends to zero, which leads to solving the following PDE

$$ \partial_{\tau} V = \frac{1}{2} \sigma^2(S)^2 S^2 \partial_{SS} V + \inf_{\alpha} \{-\alpha \partial_x V\}. $$

Finally, we give a brief outline of the numerical method used to solve the coupled PDEs (4.33)-(4.35) along with the corresponding initial and boundary conditions. We follow [75] to provide an informal discretization of the latter using a semi-Lagrangian approach. We refer the reader to the reference [49] for more details concerning the semi-Lagrangian method for HJB equations.

Along the trajectory $x = x(\tau)$ defined by

$$ dx^{\alpha} = \alpha d\tau \quad (4.37) $$

equation (4.35) can be written as

$$ \inf_{\alpha \le 0} \left\{ \frac{DV}{D\tau}(\alpha) - \mathcal{L}V - g(t, S, x, -\alpha) \right\} = 0, $$

where the operator $\mathcal{L}V$ is given by

$$ \mathcal{L}V = \frac{1}{2} \sigma^2(S) S^2 \partial_{SS} V, $$

and where the Lagrangian derivative $\frac{DV}{D\tau}(\alpha) = \partial_{\tau}V + \alpha \partial_x V$ The Lagrangian derivative is the rate of change of $V$ along the trajectory (4.37). At the same time we keep in mind that $\tilde{P}$ is the solution of the PDE (4.33) where we replace $\dot{x}$ by $-\alpha$ in the enlarged volatility.
Define a set of nodes [$s_0, s_1, ..., s_{i_{max}}$], [$x_0, x_1, ..., x_{j_{max}}$], discrete times $\tau^n = n\Delta\tau$, and localize the control candidates to values in finite interval [$\alpha_{min}, \alpha_{max}$]. Let $V(\tau^n, s_i, x_j)$ denote the exact solution to Equation (4.35), $\tilde{P}(\tau, s_i, x_j, \alpha)$ is the solution to Equation (4.33) when the control value is $\alpha \in [\alpha_{min}, \alpha_{max}]$ and $\tilde{\Delta}(\tau, s_i, x_j, \alpha) = \partial_S \tilde{P}(\tau, s_i, x_j, \alpha)$ its partial derivative w.r.t the asset price. Let $V_{i,j}^n, \tilde{P}_{i,j}^n(\alpha), \tilde{\Delta}_{i,j}^n(\alpha)$ and $\mathcal{L}_u V_j^n$ denote, respectively, the discrete approximation to the exact solution, price, its first derivative and the differential operator as in (4.32).

Let $\alpha_{i,j}^n$ denote the approximate value of the control variable $\alpha$ at mesh node $(\tau^n, s_i, x_j)$. Then the approximate $\frac{DV}{D\tau}(\alpha)$ at $(\tau^{n+1}, s_i, x_j)$ by the following:

$$ \left(\frac{DV}{D\tau}(\alpha)\right)_{i,j}^{n+1} \approx \frac{1}{\Delta\tau}(V_{i,j}^{n+1} - V_{i,\hat{j}}^n) \quad (4.38) $$

where $V_{i,\hat{j}}^n$ is an approximation of $V(\tau^n, s_i^n, x_{\hat{j}}^n)$ obtained by linear interpolation of the discrete values $V_{i,j}^n$, with $(s_i^n, x_{\hat{j}}^n)$ given by solving Equation (4.37) backwards in time for fixed $x_{i,j}^n$ to give

$$ x_{\hat{j}}^n = x_j - \alpha_{i,j}^{n+1} \Delta\tau. $$

Our final discretization is then

$$ V_{i,j}^{n+1} = \Delta\tau(\mathcal{L}_h V)_{i,j}^{n+1} + \inf_{\alpha_{i,j}^{n+1} \in [\alpha_{\min}, \alpha_{\max}]} \left\{ V_{i,j}^n + \Delta\tau \left( -\alpha_{i,j}^{n+1} \tilde{P}_{i,j}^{n+1}(\alpha_{i,j}^{n+1}) + \lambda \sigma_0 (x_j)^2 S_j^{\beta} (\tilde{\Delta}_{i,j}^{n+1}(\alpha_{i,j}^{n+1}))^2 \right) \right\}. \quad (4.39) $$

We need to solve a local optimization problem at each node at each time step in Equation (4.39). In fact, we are seeking the global minimum of the local optimization problem. If the set of controls $[\alpha_{min}, \alpha_{max}]$ is discretized with spacing $h$, then a linear search of the control space will converge to the viscosity solution of the HJB Equation (4.35) as argued in [154]. However, the uniqueness of the solution is not guaranteed with such method.

Results are presented in Figures 4.9 and 4.10 for the set of parameters given in Table 4.1 and a CEV parameter $\beta = 0.1$. The first figure gives the surface of the control variable by fixing one of the variables as a function of the remaining two others. The expected cost case does not seem to be sensitive to the moneyness as confirmed by the Heat Map in 4.10. The mean-variance case, on the contrary incites the agent to acquire more at the beginning of the strategy.

**Remarks 7.** (i) The numerical experiment is very sensitive to the option pricing. To lead the optimization, one needs both the option price and delta over each time and space grid and for each potential control $\alpha$. Both are computed using finite differences methods, which turns out to be quite slow. (ii) The option price variation over small time periods is very small. This leads to some discrete pattern of the control. (iii) More sophisticated methods need to be developed for solving such an optimization problem where the controlled function $f$ is the solution of a PDE. In particular, forward-backward BSDEs can be very useful. (iv) Finally, we can conclude from the expected cost case that, even though the impact is nonlinear in the option price, the strategy seems to be independent of the asset price evolution. This, however, needs to be invested further. The risk criterion (here the variance) is what makes the agent sensitive to the market uncertainty.
<figure><img src="image_55.png" /><figcaption>FIGURE 4.9 – The surface of the control variable $\alpha$ by fixing one of the directions (time $t$, asset price S and inventory $x$) as a function of the two remaining ones : top for $\lambda = 0$, bottom for $\lambda = 100$.</figcaption></figure>
<figure><img src="image_56.png" /><figcaption>FIGURE 4.10 – Heat maps showing the density of the trading speed and inventory : top for $\lambda = 0$, bottom for $\lambda = 100$.</figcaption></figure>
# Chapitre 5

## Range-based proxies and rough volatility

**Abstract**— In [89], it has been shown that volatility exhibits a fractional behavior with a Hurst exponent $H < 0.5$, changing the typical perception of volatility. In their study, Gatheral and co-authors used the realized volatility. In our analysis, we explore range-based proxies of the volatility process to confirm their findings on more available data (range-based) and non-standard assets. We find that the log-volatility based on range-based estimators behaves like a fractional Brownian motion with $H$ lower than 0.1. We also find that rough fractional stochastic volatility model (RFSV) is a relevant volatility model. Moreover, the prediction power of this model outperforms that of the AR, HAR and GARCH models in most cases.

**Keywords** : Range-based volatility; Garman-Klass; Parkinson; volatility scaling; fractional Brownian motion; fractional Ornstein-Uhlenbeck; volatility forecasting.

### 5.1 Introduction

Volatility plays a crucial role in many areas of finance and economics such as risk management and portfolio selection. It is known to be both time-varying and predictable, and stochastic volatility models are one way to deal with these features. As a consequence, its modeling and forecasting spurs the interest of many authors, academics and practitioners alike.

In the financial markets, the common practice is to represent asset prices by a continuous semi-martingale. A given log-price $\log S_t$ is defined by

$$d \log S_t = \mu_t dt + \sigma_t dW_t$$

where $\mu_t$ is a drift term, $W_t$ a Brownian motion and $\sigma_t$, the key ingredient is the process volatility. The Black-Scholes framework assumes the volatility to be either constant or deterministic. Such specification proved to be inadequate in the late eighties. The main reason is the inconsistency of the Black-Scholes model with the observed European options. This gave rise to alternatives such as local volatility models, such as Dupire's, see [69], and Derman and Kani's, see [63]. These models consider $\sigma_t$ as a deterministic function of time and asset price. Even though they enable us to perfectly fit a given implied surface, its dynamic is quite unrealistic. An other alternative is to model the volatility $\sigma_t$ by a continuous Brownian semi-martingale, typically correlated with $W_t$. These so-called stochastic volatility models have been the center of interest of many authors. We cite amongst such stochastic volatility models, the Hull and White model [101], the Heston model [97] and the SABR model [95]. However, generated option prices are still not consistent with observed European option prices. The reader can refer to [87] for a review of different approaches. More recently market practice is to use so-called local-stochastic volatility models which both fit the market exactly and generate reasonable
dynamics.

Since the volatility is a latent variable, the first issue one faces when trying to exhibit its statistical properties is its estimation. One can only estimate it using the underlying asset prices or quoted options. For example, when we have daily stock returns, the squared variance is a well known volatility proxy, also known as realized volatility. It measures changes on the asset return over a specified period of time. If high frequency data is available (the whole price process during the day), this proxy is more precise and is used as the daily volatility.

Access to high frequency data is sometimes costly and/or unavailable for certain assets. Therefore, other proxies are used to estimate daily volatility. If we only have closing prices and need to estimate volatility on a daily basis, we can use the squared daily returns. A compromise can be found as, in addition to closing prices, open, high and low daily prices are available for most financial data sets. In [139], Parkinson was first to introduce an advanced volatility estimator using these so-called range prices instead of just closing prices, and this enables us to overcome the issues of the first two approaches, and present a way to better estimate daily volatility.

Statistical properties of volatility estimators raise interesting questions, particularly in relation to the smoothness of the volatility process. Researchers aim to uncover the underlying mechanisms that generate the data, using the empirical scaling evidence as a stylized fact that any theoretical model should also reproduce. For example, it is common belief that volatility exhibits what is commonly known as long range dependence. The implication of this it that volatility shocks today will influence its expectation in the same direction, see [33, 72] among others.

Stochastic or local volatility models mentioned earlier assume that the smoothness of the sample path of the volatility is that of a Brownian motion (1/2 - ε Hölder continuous for any ε > 0). In [56], Comte and Renault choose to address the question of long range dependence in terms of the regularity of the driving process. Their idea was to exploit the fractional Brownian motion. Recall that a fractional Brownian motion W^H with Hurst parameter H ∈ (0, 1) is an a.s. continuous, centered, self-similar Gaussian process with stationary increments and a covariance satisfying:

$$Cov(W_t^H, W_s^H) = \frac{1}{2}(|t|^{2H} + |s|^{2H} - |t-s|^{2H}), t, s \in \mathbb{R}$$

They proposed the fractional Ornstein-Uhlenbeck volatility model with a Hurst parameter greater than 1/2 named fractional volatility model (FSV). Such model is H-ε Hölder continuous with H > 1/2. More interest grew from this model and others develop deeper analysis and calibration, see [51, 55] among others. Later, Gatheral and co-authors, see [89], challenged the previous results and established that the log-volatility process is very close to that of a fractional Brownian motion with Hurst parameter around 0.1 (< 1/2). They also developed the rough fractional stochastic volatility model (RFSV) which operates with different parameter properties than the FSV, and justified that their model better respects the volatility smile and the data properties.

In this paper, we conduct a similar study to [89] using range-based estimators to find the best forecasting model. We replicate their analysis step by step in order to revisit their finding on less standard assets. We actually find that the Hurst parameter for range-based proxies on our set of data is even lower than 0.1 and sometimes even close to 0, confirming that rough volatility hypothesis can not be refuted at this point, while further analyses allow to dismiss the hypothesis that it is generated from the FSV model (H > 1/2).
Our paper is organized as follows. In Sections 5.2 we give an overview of range-based volatility estimation before choosing one volatility proxy to work with. We conduct our statistical study in Section 5.3 where we find that log-volatility from range-based proxies behave like a fractional Brownian motion with Hurst exponent lower than 0.1. We validate the rough fractional stochastic volatility model (RFSV) introduced in [89] on our data in Section 5.4, compare its prediction power with other common models in Section 5.5, and finally conclude in Section 5.6.

## 5.2 Overview on range-based volatility estimation

In this section we review a few range-based volatility estimators and compare them to realized volatility.

We assume that the asset price over a one day period of time, $S_u$, follows a geometric Brownian motion:

$$dS_u = \mu S_u du + \sigma S_u dW_u,$$

where $u \in (t-1, t]$ is the time index between two consecutive days $t-1$ and $t$, $\mu$ is the drift, $\sigma$ the volatility considered constant along one day, and $W_u$ a standard Brownian motion. By Ito's lemma the log price $\log(S_u)$ follows a Brownian motion with drift $\mu^* = \mu - \frac{\sigma^2}{2}$ and volatility $\sigma$. During a day, it is common practice to assume that the drift is equal to zero (i.e. $\mu^* \approx 0$).

The volatility being a latent variable, one needs to estimate it at each period of time with the available asset prices. We denote by H, L, O, C the high, low, open and close prices respectively. The log-returns $r_t$ are defined by $r_t = \log(C_t) - \log(O_t)$ (or $\log(C_t) - \log(C_{t-1})$ when taking only close-to-close prices). Volatility changes, and our first interest is to be able to estimate it in a precise way and on a daily basis. Taking the squared return $r_t^2$ is one possible solution. It is an unbiased estimator of $\sigma^2$ under the normal log-returns assumption with zero mean. However, this estimator is quite noisy.

Range prices bring more consistency and information about the entire process than the close-to-close prices. In this context, Beckers, see [21], shows that volatility estimators can be improved by incorporating high and low prices, along with closing prices. Of course, range-based volatility estimators are not as efficient as realized volatility under ideal conditions, i.e. estimated from high frequency data, but remain a good alternative when this data is not available. It was shown in [4] and [151] that these range-based estimators are robust to microstructure noise and prove to be efficient and simple to compute.

In [139], Parkinson was first to develop a classical range-based estimator using high and low prices information expressed by the following formula:

$$\sigma_{\text{Parkinson}}^2 = \frac{1}{4\log2} \left( \log \frac{H_t}{L_t} \right)^2.$$

The correctness of this formula relies on a constant volatility assumption during each one day time period. It exploits the extreme value method, and the coefficient $\frac{1}{4\log2}$ is nothing but the variance of the range variable (i.e the different between the minimum and maximum). The reader can refer to [139] for details and to [73] for the asymptotic distribution of the range.

The Parkinson estimator is asymptotically unbiased under the assumption that a geometric Brownian motion without drift can describe the path of the asset price changes.
Later, Garman and Klass, see [86], established a more efficient estimator that takes the following form :

$$
\sigma_{\text{GK}}^2 = 0.511 \left( \log \frac{H_t}{L_t} \right)^2 - 0.019 \left( \log \frac{C_t}{O_t} \left( \log \frac{H_t}{O_t} - \log \frac{L_t}{O_t} \right) \right. \\
\left. - 2 \log \frac{H_t}{O_t} \log \frac{L_t}{O_t} \right) - 0.383 \left( \log \frac{C_t}{O_t} \right)^2.
$$

This estimator combines the squared return and Parkinson volatility estimators into a new estimator with smaller variance. Garman and Klass actually proved that this estimator is optimal in a mean-variance sense among a certain class of estimator, see [86].

A more practical estimator is recommended with nearly the same efficiency but eliminates the small cross-product terms expressed as :

$$
\sigma_{\text{GK}}^2 = \frac{1}{2} \left( \log \frac{H_t}{L_t} \right)^2 - (2\log 2 - 1) \left( \log \frac{C_t}{O_t} \right)^2 . \qquad (5.1)
$$

We rely on this estimator when using the Garman-Klass volatility proxy.

Because in [86] and [139] log-prices are assumed to follow geometric Brownian motion with no drift, many authors tried to correct this mismatch for securities with non-zero mean. In [145], more sophisticated drif-independent measures of volatility are introduced. The Rogers-Satchell estimator for example takes the form :

$$
\hat{\sigma}_{\text{RS}}^2 = \log \frac{H_t}{O_t} \left( \log \frac{H_t}{O_t} - \log \frac{C_t}{O_t} \right) + \log \frac{L_t}{O_t} \left( \log \frac{L_t}{O_t} - \log \frac{C_t}{O_t} \right).
$$

Kunimoto [113] and Yang-Zhang [155] also deserve to be mentioned. Yang-Zhang estimator, however, can only be used over multiple days and therefore won't be interesting for our analysis.

The previously mentioned estimators are unbiased estimators of σ². When applying the square root and estimating the volatility, all σ estimators are biased. This was expected since E[σ²] and E[σ]² are generally different.

In term of efficiency, all previous estimators exhibit very substantial improvements compared to the close-to-close estimator. Efficiency measure of a volatility estimator $\hat{\sigma}_i^2$ is defined as the ratio of the variance of this estimator and the variance of the close-to-close estimator $\hat{\sigma}_{CC}^2$ :

$$
\text{Eff}(\hat{\sigma}_t^2) = \frac{\text{Var}(\hat{\sigma}_{CC}^2)}{\text{Var}(\hat{\sigma}_i^2)}.
$$

By definition, the squared return estimator has efficiency 1. Parkinson reported that his estimator is 2.5 to 5 times more efficient than simple close-to-close variance estimator. Garman-Klass reports 7.4 while Rogers-Satchell efficiency is 6.0 and Kunimoto is 10.

In light of [13], we compare the performance and distributional properties of different range-based volatility estimators on the S&P 500 over a 3786 period from January 2000 to April 2015. Our benchmark is the realized volatility from the Oxford-Man Institute of Quantitative Finance Realized Library.
To perform the comparison we use the following measures :

The mean squared error defined by :

$$
\text{MSE}(\sigma_{\text{estimated}}) = \mathbb{E} \left[ (\sigma_{\text{estimated}} - \sigma_{\text{benchmark}})^2 \right].
$$
The mean absolute bias given by :

$$MAD(\sigma_{\text{estimated}}) = \mathbb{E}[|\sigma_{\text{estimated}} - \sigma_{\text{benchmark}}|].$$

The proportional bias expressed as :

$$\text{Prop.Bias}(\sigma_{\text{estimated}}) = \mathbb{E}\left[\left(\frac{\sigma_{\text{estimated}}}{\sigma_{\text{benchmark}}} - 1\right)\right].$$

Recall that the benchmark volatility $\sigma_{\text{benchmark}}$ here corresponds to realized volatility.

The results are given in Figures 5.1 and 5.2 and Table 5.1.

<figure><img src="image_57.png" /><figcaption>FIGURE 5.1 – Comparison between the GK proxy and the RV proxy as a benchmark. Top is the evolution of the volatility time series for both estimators. Bottom graph is a scatter plot of the GK estimator to the RV estimator.</figcaption></figure>

<table><caption>TABLEAU 5.1 – Comparison measures for different volatility estimators</caption><thead><tr><th></th><th>MSE</th><th>MAD</th><th>Prop.Bias</th><th>Std.Dev</th></tr></thead><tbody><tr><td>σ̂<sub>CC</sub></td><td>0.533 × 10<sup>-4</sup></td><td>0.0051</td><td>-0.1275</td><td>0.0092</td></tr><tr><td>σ̂<sub>P</sub></td><td>0.092 × 10<sup>-4</sup></td><td>0.0021</td><td>-0.0849</td><td>0.0063</td></tr><tr><td>σ̂<sub>GK</sub></td><td>0.094 × 10<sup>-4</sup></td><td>0.0021</td><td>-0.1313</td><td>0.0058</td></tr><tr><td>σ̂<sub>RS</sub></td><td>0.206 × 10<sup>-4</sup></td><td>0.0028</td><td>-0.1762</td><td>0.0062</td></tr></tbody></table>

We can see in Table 5.1, Figures 5.1 and 5.2 the following :

- Range-based estimators are lower than the benchmark.

- Range-based estimators reduce the variance compared to squared returns.

- The Garman-Klass estimator is the closest to intraday realized variance and has the smallest variance.
<figure><img src="image_58.png" /><figcaption>FIGURE 5.2 – Comparison between different range-based estimators.</figcaption></figure>

— Range-based estimators seem to exhibit the same regularity as that of realized volatility.

We would like to confirm the last point. In our study, we focus on the volatility based on the Garman-Klass proxy. Our goal is to confirm that the scaling properties discussed in [89], are also satisfied by the volatility when using these proxies.

## 5.3 Range-based volatility as spot volatility proxy : empirical results

### 5.3.1 The scaling of the Garman-Klass proxy

We carry out our analysis on the volatility proxy for a variety of assets. In [89], Gather and co-authors use common indexes for which high frequency based realized volatility is available on the Oxford-Man Institute of Quantitative Finance Realized Library (S&P 500, Bund ...). In this paper, we choose to apply our analysis on more "exotic" assets (S&P 400, IBEX 35, IBOV, S&P 100, INDU, SHSZ300, MEXBOL, FTSE 100, XIN9I, HSI), and some stocks (TOTAL, ASX200, GOOGLE and MICROSOFT). Most of these assets are not available on the Library, and more importantly, since many financial institutions still do not have access to high frequency based proxies, the choice of range-based is very convenient. Range data availability allows for their computation for any class of assets and any assets tickers.

We choose to present the analysis for the S&P 100¹ and IBEX 35², and give numerical results for the

<footnote n="1">1. The S&P 100 is a subset of the S&P 500 and includes 102 leading U.S. stocks with exchange-listed options. Constituents of the S&P 100 represent about 63% of the market capitalization of the S&P 500 and almost 51% of the market capitalization of the U.S. equity markets as of January 2017.</footnote>

<footnote n="2">2. IBEX 35 is the benchmark stock market index of the Bolsa de Madrid, Spain's principal stock exchange. It is a market capitalization weighted index comprising the 35 most liquid Spanish stocks traded in the Madrid Stock Exchange General Index</footnote>
remaining ones. We focus mainly on the Garman-Klass proxy and present results on the Parkinson estimator for verification.

The set of data corresponds to 2521 trading days from April 19, 2005 to April 22, 2015. Volatility proxies are based on range-data extracted from Bloomberg database. Let $\sigma_{t_0}, \sigma_{t_1}, ..., \sigma_{t_N}$, be the time series of the GK proxy computed for this period, where $t_{i+1} - t_i$ corresponds to one business day.

Our scaling measure $m(q, \Delta)$ is the q-th absolute moment of the increments of log-volatility and are defined by:

$$m(q, \Delta) := \frac{1}{N} \sum_{k=1}^{[N/\Delta]} |\log(\sigma_k \Delta) - \log(\sigma_{(k-1)\Delta})|^q,$$

for different $q > 0$ and lags $\Delta$ going from 1 to about 400 days. Our goal is to revisit the finding in [89] that the spot log-volatility has the same scaling properties as a fractional Brownian motion with Hurst exponent $H < 1/2$, and therefore one can model it with such process. It is worthy to mention to following remarks:

**Remark 11.** — The quantity $m(q, \Delta)$ is the discrete equivalent of $\mathbb{E}[|\log(\sigma_{\Delta}) - \log(\sigma_0)|^q]$. Recall that the fractional Brownian motion $W^H$ verifies:

$$\mathbb{E}[|W_{t+\Delta}^H - W_t^H|^q] = \tilde{K}_q \Delta^{qH},$$

— We are expecting the volatility to behave closely to the fractional Brownian motion. As a result, we would observe the following relationship:

$$m(q, \Delta) \sim K_q \Delta^{qH}. \qquad (5.2)$$

— Given that data is finite and not time-equidistant (unavailable data on weekends for example), the measure $m(q, \Delta)$ is the result of averaging over all possible increments by taking a rolling window, and selecting only increments that correspond to the chosen lag $\Delta$ (log-vol increments between two successive volatility measures between Friday through Monday are considered as 3 days lag).

To verify the validity of (5.2), we plot $\log(m(q, \Delta))$ against $\log(\Delta)$ for different values of $q$. Depending on the results of this first regression, we can write, for a given $q$,

$$m(q, \Delta) \sim K_q \Delta^{\zeta_q}, \qquad (5.3)$$

where $\zeta_q$ defines a general scaling function.

Results are displayed in Figures 5.3 and 5.4. We can notice the following:

— The values of $m(q, \Delta)$ for different $q$ against $\log(\Delta)$ lie within a straight line (left figures). This confirms that both S&P 100 and IBEX 35 exhibit a scaling property given by Equation (5.3).

— The R-squared values given in Table 5.2 confirm that the data is close to the fitted regression line for all values of $q$.

— Plotting $\zeta_q$ as a function of $q$ (right), confirms our expectation; the scaling is linear in $q$ and verifies Equation (5.2), with $H = 0.081$ for S&P 100 and $H = 0.072$ for IBEX 35.
<figure><img src="image_59.png" /><figcaption>FIGURE 5.3 – log $m(q, \Delta)$ as a function of log $\Delta$ (left), $\zeta_q$ (blue) and $0.084 \times q$ (green) (right), S&P 100 (Garman Klass volatility).</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>FIGURE 5.4 – log $m(q, \Delta)$ as a function of log $\Delta$ (left), $\zeta_q$ (blue) and $0.072 \times q$ (green) (right), IBEX35 (Garman Klass volatility).</figcaption></figure>

<table><caption>TABLEAU 5.2 – The R-squared measure of the regression $log(m(q, \Delta)) \sim log(\Delta)$</caption><thead><tr><th>q</th><th>0.5</th><th>1</th><th>1.5</th><th>2</th><th>3</th></tr></thead><tbody><tr><td>S&P 100</td><td>0.9640</td><td>0.9760</td><td>0.9800</td><td>0.9797</td><td>0.9717</td></tr><tr><td>IBEX 35</td><td>0.9225</td><td>0.9290</td><td>0.9287</td><td>0.9252</td><td>0.9111</td></tr></tbody></table>

Of course the $\zeta_q$ and Hq are not perfectly matched. One possible reason for this mismatch is using discrete samples. In fact, simulating fBm using the same number of points results in a slight concave figure, and ends up in recovering a Hurst parameter that slightly overestimates the real one.

In order to make sure our estimations of H do not depend on the time interval, we split the data into two periods with the same lag and re-estimate H for each period separately. The aim of this experience is to confirm that the scaling is time independent for all assets. Regressing $log(m(q, \Delta))$ on $log(\Delta)$ for each ticker and for $q = 0.5, 1, 1.5, 2, 3$ and $\Delta = 1, ..., 410$, we find that $\zeta_q$ is linear on $q$ for the GK proxy. As we can see in Table 5.6, the Hurst parameter lies between 0.01 and 0.082 confirming that the volatility process is rough. One might however notice that splitting the data resulted in a first half with a slightly greater H than the second half. We think that it is due to the presence of the 2008 crisis in the first half.
<table><caption>TABLEAU 5.3 – Estimates of H on the whole period and over two different time intervals for different indexes and stocks (Garman Klass volatility)</caption><thead><tr><th>Ticker</th><th>H for the whole period</th><th>H (first half)</th><th>H (second half)</th></tr></thead><tbody><tr><td>SP100</td><td>0.0841</td><td>0.0897</td><td>0.0714</td></tr><tr><td>IBEX35</td><td>0.072</td><td>0.0753</td><td>0.071</td></tr><tr><td>HSI</td><td>0.0516</td><td>0.0605</td><td>0.0394</td></tr><tr><td>MEXBOL</td><td>0.0627</td><td>0.0737</td><td>0.0463</td></tr><tr><td>FTSE100</td><td>0.0751</td><td>0.0708</td><td>0.0728</td></tr><tr><td>ASX200</td><td>0.0489</td><td>0.0476</td><td>0.0415</td></tr><tr><td>TOTAL</td><td>0.0774</td><td>0.0835</td><td>0.0687</td></tr><tr><td>XIN9I</td><td>0.0674</td><td>0.0649</td><td>0.069</td></tr><tr><td>SHSZ300</td><td>0.0689</td><td>0.0718</td><td>0.0636</td></tr><tr><td>BCOM</td><td>0.014</td><td>0.0099</td><td>0.0238</td></tr><tr><td>INDU</td><td>0.0804</td><td>0.0838</td><td>0.067</td></tr><tr><td>USDEUR</td><td>0.0353</td><td>0.0393</td><td>0.0321</td></tr><tr><td>IBOV</td><td>0.0685</td><td>0.0724</td><td>0.0609</td></tr><tr><td>MICROSOFT</td><td>0.06</td><td>0.0717</td><td>0.0401</td></tr><tr><td>GOOGLE</td><td>0.0656</td><td>0.0724</td><td>0.0542</td></tr><tr><td>SP400</td><td>0.0715</td><td>0.0753</td><td>0.0592</td></tr></tbody></table>

### 5.3.2 The scaling of the Parkinson volatility proxy

To ensure that the results apply to other range-based estimators, we reproduce the same analysis to the Parkinson volatility proxy, based on the same data and on the same period. Results are expressed in Figures 5.5 and 5.6.

<figure><img src="image_61.png" /><figcaption>FIGURE 5.5 – $\log m(q, \Delta)$ as a function of $\log \Delta$ (left), $\zeta_q$ (blue) and $0.082 \times q$ (green) (right), S&P 100 (Parkinson volatility).</figcaption></figure>

The scaling is again linear in $q$ with $H = 0.082$ for the S&P 100 and $H = 0.064$ for the IBEX 35. Values are very close to those found within the GK proxy. This confirms again that volatility based on the Parkinson proxy is rough.

The smoothness parameter H is detailed in Table 5.4 for all the assets for the Parkinson volatility. One more time, we give the value for the whole period, split the period into two halves and compute H for each half.
<figure><img src="image_62.png" /><figcaption>FIGURE 5.6 – log *m*(*q*, Δ) as a function of log Δ (left), ζ<sub>*q*</sub> (blue) and 0.064 × *q* (green) (right), IBEX 35 (Parkinson volatility).</figcaption></figure>

The same conclusions apply to the Parkinson volatility :

<table><caption>TABLEAU 5.4 – Estimates of H on the whole period and over two different time intervals for different indexes and stocks (Parkinson volatility)</caption><thead><tr><th>Ticker</th><th>H for the whole period</th><th>H (first half)</th><th>H (second half)</th></tr></thead><tbody><tr><td>SP100</td><td>0.0822</td><td>0.0888</td><td>0.0737</td></tr><tr><td>IBEX35</td><td>0.0644</td><td>0.0682</td><td>0.0648</td></tr><tr><td>HSI</td><td>0.0452</td><td>0.0555</td><td>0.0336</td></tr><tr><td>MEXBOL</td><td>0.0638</td><td>0.0738</td><td>0.0489</td></tr><tr><td>FTSE100</td><td>0.0774</td><td>0.0823</td><td>0.0669</td></tr><tr><td>ASX200</td><td>0.0513</td><td>0.0511</td><td>0.0422</td></tr><tr><td>TOTAL</td><td>0.0738</td><td>0.0856</td><td>0.0608</td></tr><tr><td>XIN9I</td><td>0.0595</td><td>0.0592</td><td>0.0593</td></tr><tr><td>SHSZ300</td><td>0.0591</td><td>0.0668</td><td>0.0499</td></tr><tr><td>BCOM</td><td>0.0127</td><td>0.00623</td><td>0.0251</td></tr><tr><td>INDU</td><td>0.08</td><td>0.08</td><td>0.0707</td></tr><tr><td>USDEUR</td><td>0.0265</td><td>0.0276</td><td>0.0261</td></tr><tr><td>IBOV</td><td>0.0694</td><td>0.0746</td><td>0.0603</td></tr><tr><td>MICROSOFT</td><td>0.0584</td><td>0.0685</td><td>0.0414</td></tr><tr><td>GOOGLE</td><td>0.0603</td><td>0.063</td><td>0.0542</td></tr><tr><td>SP400</td><td>0.0757</td><td>0.0822</td><td>0.0623</td></tr></tbody></table>

— H remains between 0.01 and 0.09 with most assets around 0.07. This fact confirms that volatility is rough.

— H is higher for the first period. This might be explained by the fact that this period (period 2005-2010) contains the 2008 crisis.

— H for Parkinson proxy of the BCOM<fnref n="3" /> asset is almost 0.

<footnote n="3">3. BCOM corresponds to the B Communications Ltd, which is a publicly traded holding company, headquartered in Israel, whose sole asset is a controlling interest in Israeli telecommunications provider Bezeq.</footnote>
### 5.3.3 Distribution of the increments of the log-volatility

Now that we have established the common scaling behavior for our volatility proxies on the given data, we will focus on the S&P 100 for the following results. Certainly, we will ensure that the same results are common for other assets, but unless specified otherwise, all the plots concern the S&P 100.

It is well-known that the increments of log-volatility distribution is very close to the normal distribution, see for example [10]. This is also what we find in our data, see Figure 5.7. Moreover, rescaling the density by $\Delta^H$ for any given lag $\Delta$ recovers the 1-day increments density. This is consistent with the fractional Brownian motion with Hurst parameter H as seen in the previous section.

<figure><img src="image_63.png" /><figcaption>FIGURE 5.7 – Histograms for various lags $\Delta$ of the increments $\log(\sigma_{t+\Delta}) - \log(\sigma_t)$ of the Garman-Klass S&P 100 log-volatility; the normal fit to distribution of the $\Delta$-days increments (red); normal fit to the 1-day increments rescaled by $\Delta^H$ (dashed blue)</figcaption></figure>

## 5.4 RFSV model validation using range-based proxies

In the following section, we test the accuracy of the Rough fractional volatility model (RFSV) introduced in [89] using the volatility data based on the Garman-Klass proxy. We will ensure that the model reproduces the same behavior as the data.
### 5.4.1 The model

Empirical results of Section 5.2, show that the increments of log-volatility based on range proxies for various assets appear to have a scaling property with constant Hurst parameter. We also made sure, through Section 5.3.3 that their distribution is close to a normal distribution. Based on these results, the log-volatility increments can be modeled by the increments of a fBm as the following:

$$ \log \sigma_{t+\Delta} - \log \sigma_t = v(W_t^H - W_{t+\Delta}^H), \quad (5.4) $$

where $W^H$ is a fractional Brownian motion with the Hurst parameter estimated through the scaling of the volatility and $v$ is a positive constant corresponding to the volatility of the increments. We can rewrite Equation (5.4) under the form:

$$ \sigma_t = \sigma \exp(vW_t^H), $$

where $\sigma$ is a positive constant.

One of the drawbacks of this model is that it is not stationary. As a matter of fact, stationarity is a property that is desirable and useful for modeling time series. A possible model that keeps this property along with the fractional scaling is the fractional Ornstein-Uhlenbeck (fOU in short) process with a very long mean-reversion.

The fractional Ornstein-Uhlenbeck process $X_t$ is a stochastic process satisfying the stochastic differential equation:

$$ dX_t = -\alpha X_t dt + \nu dZ_t, \quad X_0 = 0 $$

where both $\nu$ and $\alpha$ are positive parameters. When $Z_t$ is the standard Brownian motion, we get the standard Ornstein-Uhlenbeck, see [136]. Our interest, however, is the when $Z_t = W_t^H$. We also consider an arbitrary initial point $X_0 = m$ instead of 0. The SDE followed by the process of log-volatility becomes:

$$ dX_t = -\alpha(X_t - m)dt + \nu dW_t^H, \quad x_0 = m \quad (5.5) $$

where $m \in \mathbb{R}$ and $(W_t^H)$ is the fBm with Hurst parameter $H$.

Equation (5.5) is then solved using the following explicit representation:

$$ X_t = v \int_{-\infty}^{t} e^{-\alpha(t-s)} dW_t^H + m \quad (5.6) $$

where the stochastic integral with respect to fBM is simply a path-wise Riemann-Stieljes integral ([51]). Lastly, we recover the volatility, and thus define the RFSV model on the time interval $[0,T]$:

$$ \sigma_t = \exp(X_t), \quad t \in [0, T], $$

where $(X_t)$ satisfies equation (5.6) for some $\nu > 0$, $\alpha > 0$, $m \in \mathbb{R}$ and $H < 1/2$ the measured smoothness of the volatility. In addition to the stationarity of such a model, choosing $\alpha \ll 1/T$ allows the log-volatility to behave locally (at time scales smaller than $T$) as a fBm. This observation is formalized by Proposition 3.1 in [89] we recall below:

**Proposition 5.** Let $W^H$ be a fBm and $X^\alpha$ defined by (5.6) for a given $\alpha > 0$. As $\alpha$ tends to zero,

$$ \mathbb{E}[\sup_{t \in [0,T]} |X_t^\alpha - X_0^\alpha - \nu W_t^H|] \to 0 $$
**Remark 12.**:

- Proposition (5) implies that within the interval $[0,T]$, and under the condition $\alpha \ll 1/T$, we can proceed as if the log-volatility process were a fBm. Setting $\alpha = 0$ allows to recover the simple non-stationary fBm (5.4).

- The RFSV differs from the classical FSV model of Comte and Renault, see [56], in that, instead of taking $H > 1/2$ and $\alpha$ large in FSVr, the RFSV model is defined for $H < 1/2$ and $\alpha$ small (actually $\alpha$ is chosen not to be equal to 0 only so that the volatility satisfies stationarity).

- The choice $H < 1/2$ is consistent with both the statistical properties of the data and generates a term structure of volatility skew that matches the observations.

- The choice of the fOU is for convenience and simplicity. Other models that imitate the fBm behavior at reasonable time scales and are stationary can be considered as well.

- The RFSV process reproduces approximately the exact scaling property as the fBm. This is a consequence of the following corollary:

**Corollary 1.** Let $q > 0$, $t > 0$, $\Delta > 0$. As $\alpha$ tends to zero, we have:

$$ \mathbb{E}[|X_{t+\Delta}^{\alpha} - X_{t}^{\alpha}|^q] \to v^q K_q \Delta^{qH}. $$

A detailed comparison between the RFSV and FSV can be found in [89].

### 5.4.2 Model validation

In previous sections, we come to the conclusion that RFSV model seems to be a relevant volatility model based on empirical results. But the question that arises is whether the estimated range-based volatility proxies, from simulated data with RFSV as the volatility process, behaves like the underlying process (in terms of the scaling properties). Our goal is to investigate this question. To do so, we simulate the spot volatility process using the RFSV model, simulate intraday prices, recover range prices (open, close, high and low), and finally estimate the range volatility from the simulated range prices. We will be able, on the one hand, to compare the behavior of the known real spot volatility process to that of the proxy used for its approximation. We will also be able to estimate the realized volatility and compare it to the Garman-Klass and to real data.

Spot volatility is simulated using RFSV model for 2,521 days. Since range-based volatility assumes the volatility process to be constant within the day, we take into account the randomness of the intraday prices $P_u$ where $u$ frequency is of the order of a few seconds, and simulate the spot volatility $\sigma_t$ on a daily basis. We choose parameters that are consistent with the S&P 100 Garman-Klass volatility (even though a precise estimation for $v \alpha$ and $m$ is not simple) and consistent with our empirical estimates from Section 1.3, i.e. $H = 0.08$, $v = 0.3$, $m = X_0 = -5$ and $\alpha = 5 \times 10^{-4}$.

To simulate the volatility and price paths, we proceed as the following:

- Simulate fBm using a wavelet-based synthesis, see [2].

- Simulate the log-volatility process $X$ for each day $n$ according to a discrete scheme:

$$ X_{n+1} - X_n = v(W_{n+1}^H - W_n^H) + \alpha(m - X_n). $$

- Simulate the asset price $P$ by taking:

$$ P_{n+(j+1)\delta} - P_{n+j\delta} = P_{n+j\delta}\sigma_n\sqrt{\delta}U_j, $$

where the $U_j$ are iid standard Gaussian variables.
— Extract range prices (open, close, high and low prices) for each day, and compute the realized volatility and the Garman-Klass range-volatility.

We present in Figure 5.8 a plot of the Garman-Klass proxy from S&P 100 data along with the simulated spot volatility described above.

We compare the plots of the S&P 100 Garman-Klass proxy with the simulated path. Graphically, it seems that estimated volatility exhibits the same behavior as the simulation, at least to a visual extent. This was already verified in [89]. A zoom in or a zoom out gives typically the same kind of qualitative properties. To compare the smoothness of the real spot volatility (simulated paths), with that of the

<figure><img src="image_64.png" /><figcaption>FIGURE 5.8 – Garman Klass volatility of S&P 100 (above) and simulated paths (below)</figcaption></figure>

Garman-Klass proxy recovered from simulated prices, we repeat the analysis of Section 5.2. We plot in Figure 5.9 $log(m(q, \Delta))$ as a function of $log(\Delta)$ for both the real volatility and volatility proxy. The Hurst exponent of the estimated Garman-Klass proxy is relatively close to the true one ($H = 0.079$ compared to $H = 0.056$). The mismatch can in particular be due to simulation bias.

<figure><img src="image_65.png" /><figcaption>FIGURE 5.9 – $log(m(q, \Delta))$ as a function of $log(\Delta)$ for the real spot volatility with $H = 0.08$ and the Garman Klass proxy based on simulations.</figcaption></figure>
5.4.3 FSV vs. RFSV

We have shown from our analysis that the RFSV model hypothesis can not be rejected. In light of
empirical results on data, we would like to test if we can reject the FSV model through simulating
the volatility process. To do so, we analyze the behavior of the smoothing function *m*(*q*, *Δ*) for small
and large lags, using simulated volatility and asset prices, and see if the Garman-Klass and realized-
volatility based on these simulations, behave like the one found for real data.

The first basic difference between RFSV and FSV is the range of the Hurst exponent values for the frac-
tional Brownian motion. First, we consider the simple non mean-reverting fractional volatility model
with H > 0.5, i.e. σ_t = σ_0 e^(VWt^H) with H = 0.7 and η = 0.25. We take a look at the scaling behavior of the
realized-volatility and Garman-Klass volatility proxies based on this model in Figure 5.10. In Figure

<figure><img src="image_66.png" /><figcaption>FIGURE 5.10 – log(m(q, Δ)) as a function of log(Δ) for the data (left) Realized Volatility (center) and Garman-Klass volatility (right) over the simulated paths for q = 1, 2, 3, 4, 5</figcaption></figure>

5.10, we see that the scaling of log(Δ) is very close to a straight line for the RV and GK estimators. The resulting smoothing parameter found in these figures is close to the original one (H = 0.69 for RV estimator and H = 0.64 for GK estimator).

We are aware that such model, leading to crazy volatility values, does not make sense without mean-
reversion. However, this allows us to exclude a fractional volatility model of the form σ_t = σ_0 e^(VWt^H) with
H > 0.5.

Unlike the RSFV model where the mean-reversion is intrinsic to the model for α = 0, Comte and Re-
nault impose α to be large enough, i.e. α ≫ 1/T where T is the time horizon of interest, to verify this
property. We would like to test to what extent this model can be misleading in estimating the smooth-
ness of the diffusion process. We compare the FSV and RFSV models for the set of parameters given in
Table 5.5, which leads to the simulated time series given in Figure 5.11.
<table><caption>TABLEAU 5.5 – Parameters values used for simulating the FSV and RFSV models and the mean and variance of the simulated time series</caption><thead><tr><th></th><th>FSV</th><th>RFSV</th></tr></thead><tbody><tr><td>H</td><td>0.7</td><td>0.08</td></tr><tr><td>α</td><td>0.25</td><td>5 × 10<sup>-4</sup></td></tr><tr><td>ν</td><td>0.25</td><td>0.45</td></tr><tr><td>m</td><td>-4.5</td><td>-5</td></tr><tr><td>X<sub>0</sub></td><td>-4.5</td><td>-5</td></tr><tr><td>E[log(σ)]</td><td>-4.6</td><td>-4.7</td></tr><tr><td>Var[log(σ)]</td><td>0.21</td><td>0.33</td></tr></tbody></table>

<figure><img src="image_67.png" /><figcaption>FIGURE 5.11 – Plot of one path of the FSV and RFSV models for H = 0.7 and H = 0.08 observed daily (top two), and every 10 days (bottom two).</figcaption></figure>

In Figure 5.11, we see that both processes satisfy mean-reversion. The roughness of the RFSV model is quite clear when we observe the process daily, however, when observed every 10 days, it is less obvious to say which is the FSV and which is the RFSV.

In Figure 5.12, we check one more time the scaling behavior of the two processes and their RV and GK proxies estimated on a 24 hours windows and 1 second observed prices. We confirm through Figure 5.12 the following key results :

- The smoothing function $\zeta_q$ for the RFSV model keeps the same pattern for short and long time scales, $(\log(m(q, \Delta)))$ is close almost linear w.r.t $\log(\Delta))$.
<figure><img src="image_68.png" /><figcaption>FIGURE 5.12 – log(*m*(*q*, *Δ*)) as a function of log(*Δ*) for the simulated paths (left), Realized Volatility (center) and Garman Klass volatility (right) over the simulated paths for *q* = 1, 2, 3, 4, 5, for FSV model (top) and RFSV model (bottom)</figcaption></figure>

- The FSV model seems to exhibit two slopes. At small scales, the slope is close to that of the Hurst exponent of the fractional Brownian motion that drives the process, i.e. H ≈ 0.7. At large scales, the slope gives a value close 0. Actually, the stationarity of the process at large scales is responsible for such estimation.

Finally, in order to verify the impact of discretization on the estimators, we consider that prices are observed on an 8 hours time window every 1, 5 or 10 minutes. We compute the RV and GK estimators based on these observations. Results are shown in Figure 5.13 for the FSV (this test on the RFSV model does not bring more information).

<figure><img src="image_69.png" /><figcaption>FIGURE 5.13 – log(*m*(*q*, *Δ*)) as a function of log(*Δ*) for realized volatility (top) and Garman-Klass volatility (bottom) based on an 8 hours window with 1min (right), 5min (center) and 10min (left) discretization for the mean-reverting FSV model with H = 0.07.</figcaption></figure>
We observe an other phenomena for the GK estimator. The slope on very small scales ($\Delta = 1, 2$) gives a relatively smaller smoothing parameter (H around 0.19) than for intermediate lags (H around 0.24). We believe that this is due to the noise of this estimator. Moreover, the value of the smoothing parameter is still very small compared the one used to simulate the process.

Finally, we can conclude that FSV volatility with $H > 0.5$ behaves differently than the data. It may therefore be excluded from being a good volatility model. Rough fractional model on the other hand, seems at this point, the most plausible to model the volatility. Even though it was not illustrated in here, mean-reversion models with standard Brownian diffusion is even less plausible.

## Quantitative justification

In this section, we would like to quantify the phenomena encountered previously. Our goal is to see how estimating the smoothing parameter is affected by the lag range. We denote by $X_t$ the asset log price $X_t = \log(P_t)$. Since spot volatility does not really make sense at the intraday level (beyond seasonality), and to avoid the smoothing issue of realized volatility estimation, we assume that the volatility is constant within a day, i.e. for day $i$ and for each time $t \in [i, i+1)$ $\sigma_t = \sigma_i = \text{constant}$, with $\sigma_i = e^{\nu W_i^H}$ (although this model does not make much sense for $H > 0.5$).

Assume we have $n$ observations for each day $i$ for the price process $P_i^{j=1,...,n}$, and let $\Delta_j^n X$ the log price increments of day $i$:

$$ \Delta_j^n X = \log(P_i^j) - \log(P_i^{j-1}). $$

Using the central limit theorem (CLT) for realized volatility we have:

$$ \begin{align*} (\sigma_i^{RV})^2 &= \sum_{j=1}^{n} (\Delta_j^n X)^2 \\ &\approx \int_i^{i+1} e^{2\nu W_i^H} dt + \frac{1}{\sqrt{n}} \sqrt{2 \int_i^{i+1} e^{4\nu W_i^H} dt} \xi \quad \text{where } \xi \sim \mathcal{N}(0, 1) \\ &\approx e^{2\nu W_i^H} \left(1 + \sqrt{\frac{2}{n}} \xi\right), \end{align*} $$

which leads to the approximation:

$$ \log(\sigma_i^{RV}) \approx \nu W_i^H + \sqrt{\frac{1}{2n}} \xi \quad \text{where } \xi \sim \mathcal{N}(0, 1). $$

Taking the increments of the log volatility between $i$ and $i + \Delta$ for a given time lag $\Delta$, we have:

$$ \log(\sigma_{i+\Delta}^{RV}) - \log(\sigma_i^{RV}) \approx \underbrace{\nu(W_{i+\Delta}^H - W_i^H)}_{o(\nu\Delta^H)} + \underbrace{\sqrt{\frac{1}{2n}}(\xi + \xi')}_{o(\sqrt{\frac{1}{n}})} $$

where $\xi$ and $\xi'$ are i.i.d Gaussian variables.

This equivalence leads to the following observations:

- When the vol of vol $\nu$ is small, i.e. such $\nu^2 \Delta^{2H} \ll \frac{1}{n}$, the noise takes the upper hand. As a result, we would observe a slope close to horizontal when plotting $\log(m(q, \Delta))$ against $\log(\Delta)$.
— On the contrary, when v is large, the first term is predominant and the slope is close to qH which allows for a rather more precise estimation of the true Hurst exponent.

— In the intermediate case, two slopes can be observed; one corresponding to estimation noise, and the other proportional to the true Hurst exponent.

We illustrate these cases in Figure 5.14 below.

<figure><img src="image_70.png" /><figcaption>FIGURE 5.14 – Different cases encountered for the estimation of the slope of log(m(q, Δ)) against log(Δ).</figcaption></figure>

## 5.5 Forecasting range-based volatility using the RFSV model

The purpose of this section is to compare the predictability of the RFSV model with other commonly used models such as the AR, HAR or GARCH models.

### 5.5.1 Forecasting log-volatility

The key formula on which the prediction method is based is the following one :

$$ \mathbb{E}[W_{t+\Delta}^{H} | \mathcal{F}_{t}] = \frac{\cos(H\pi)}{\pi} \Delta^{H+1/2} \int_{-\infty}^{t} \frac{W_{s}^{H}}{(t-s+\Delta)(t-s)^{H+1/2}} ds, $$

where W^H is a fBm with H < 1/2 and F_t the filtration it generates, see Theorem 4.2 of [135]. By construction, over any reasonable time scale of interest, as formalized in Corollary 1, we may approximate the fOU volatility process in the RFSV model as log σ_t^2 ≈ 2vW_t^H + C for some constants v and C. Our prediction formula for the log-variance then follows :

$$ \mathbb{E}[\log \sigma_{t+\Delta}^2 | \mathcal{F}_t] = \frac{\cos(H\pi)}{\pi} \Delta^{H+1/2} \int_{-\infty}^{t} \frac{\log \sigma_s^2}{(t-s+\Delta)(t-s)^{H+1/2}} ds. \quad (5.7) $$

This formula, or rather its approximation through a Riemann sum (we assume in this section that volatilities are perfectly observed, although they are in fact estimated), is used to forecast the log-volatility 1,5 and 20 days ahead (Δ = 1,5,20).

In the spirit of [58], we compare the predictive power of Formula (5.7) with that of AR, HAR and GARCH forecasts. Recall that for a given integer p > 0, the AR(p) and HAR predictors take the following form (where the index i runs over the series of daily volatility estimates) :

— AR(p) :

$$ \overline{\log(\sigma_{t+\Delta}^2)} = K_0^\Delta + \sum_{i=0}^{p} C_i^\Delta \log(\sigma_{t-i}^2). $$
— HAR :

$$ \widehat{\log(\sigma_{t+\Delta}^2)} = K_0^\Delta + C_0^\Delta \log(\sigma_t^2) + C_5^\Delta \frac{1}{5} \sum_{i=0}^p C_i^\Delta \log(\sigma_{t-i}^2) + C_{20}^\Delta \frac{1}{20} \sum_{i=0}^{20} \log(\sigma_{t-i}^2). $$

We estimate AR and HAR using a rolling time window of 500 days. For the HAR case, we use standard linear regression to estimate the coefficients as explained in [58]. In the sequel, we consider $p = 5$ and $p = 10$ in the AR formula. Indeed, these parameters essentially give the best results for the horizons at which we wish to forecast the volatility (1, 5 and 20 days). For each day, we forecast volatility for five different indexes.

We then assess the quality of the various forecast by computing the ratio P between the mean squared error of our predictor and the approximated variance of the log-variance :

$$ P = \frac{\sum_{k=500}^{N-\Delta} (\log(\sigma_{k+\Delta}^2) - \widehat{\log(\sigma_{t+\Delta}^2)})^2}{\sum_{k=500}^{N-\Delta} (\log(\sigma_{k+\Delta}^2) - \mathbb{E}[\log(\sigma_{k+\Delta}^2)])^2}, $$

where $\mathbb{E}[\log(\sigma_{k+\Delta}^2)]$ denotes the empirical mean of the log-variance over the whole period.

We present in Table 5.6 the ratio P for different models in order to compare the RFSV prediction power with other autoregressive models for predicting the log variance.

<table><thead><tr><th>Ticker</th><th>AR(5)</th><th>AR(10)</th><th>HAR(3)</th><th>RFSV</th></tr></thead><tbody><tr><td>SP100 Δ = 1</td><td>0.451</td><td>0.446</td><td>0.443</td><td>0.466</td></tr><tr><td>SP100 Δ = 5</td><td>0.644</td><td>0.635</td><td>0.546</td><td>0.557</td></tr><tr><td>SP100 Δ = 21</td><td>0.897</td><td>0.894</td><td>0.734</td><td>0.718</td></tr><tr><td>IBEX35 Δ = 1</td><td>0.594</td><td>0.594</td><td>0.582</td><td>0.622</td></tr><tr><td>IBEX35 Δ = 5</td><td>0.843</td><td>0.824</td><td>0.728</td><td>0.728</td></tr><tr><td>IBEX35 Δ = 21</td><td>1.18</td><td>1.17</td><td>0.943</td><td>0.908</td></tr><tr><td>HSI Δ = 1</td><td>0.529</td><td>0.523</td><td>0.513</td><td>0.52</td></tr><tr><td>HSI Δ = 5</td><td>0.647</td><td>0.633</td><td>0.575</td><td>0.577</td></tr><tr><td>HSI Δ = 21</td><td>0.805</td><td>0.801</td><td>0.665</td><td>0.671</td></tr><tr><td>MEXBOL Δ = 1</td><td>0.572</td><td>0.567</td><td>0.553</td><td>0.589</td></tr><tr><td>MEXBOL Δ = 5</td><td>0.731</td><td>0.709</td><td>0.648</td><td>0.645</td></tr><tr><td>MEXBOL Δ = 21</td><td>0.922</td><td>0.917</td><td>0.757</td><td>0.764</td></tr><tr><td>FTSE100 Δ = 1</td><td>0.474</td><td>0.465</td><td>0.463</td><td>0.476</td></tr><tr><td>FTSE100 Δ = 5</td><td>0.627</td><td>0.614</td><td>0.545</td><td>0.545</td></tr><tr><td>FTSE100 Δ = 21</td><td>0.859</td><td>0.855</td><td>0.699</td><td>0.688</td></tr><tr><td>ASX200 Δ = 1</td><td>0.536</td><td>0.524</td><td>0.524</td><td>0.527</td></tr><tr><td>ASX200 Δ = 5</td><td>0.658</td><td>0.652</td><td>0.577</td><td>0.573</td></tr><tr><td>ASX200 Δ = 21</td><td>0.806</td><td>0.793</td><td>0.707</td><td>0.688</td></tr><tr><td>TOTAL Δ = 1</td><td>0.540</td><td>0.534</td><td>0.527</td><td>0.558</td></tr><tr><td>TOTAL Δ = 5</td><td>0.720</td><td>0.704</td><td>0.640</td><td>0.636</td></tr><tr><td>TOTAL Δ = 21</td><td>1.008</td><td>1.015</td><td>0.809</td><td>0.789</td></tr><tr><td>XIN9I Δ = 1</td><td>0.587</td><td>0.58</td><td>0.568</td><td>0.582</td></tr><tr><td>XIN9I Δ = 5</td><td>0.712</td><td>0.695</td><td>0.637</td><td>0.641</td></tr><tr><td>XIN9I Δ = 21</td><td>0.913</td><td>0.918</td><td>0.762</td><td>0.758</td></tr><tr><td>SHSZ300 Δ = 1</td><td>0.574</td><td>0.568</td><td>0.56</td><td>0.572</td></tr><tr><td>SHSZ300 Δ = 5</td><td>0.707</td><td>0.695</td><td>0.634</td><td>0.634</td></tr><tr><td>SHSZ300 Δ = 21</td><td>0.896</td><td>0.904</td><td>0.772</td><td>0.753</td></tr></tbody></table>
<table><caption>TABLEAU 5.6 – Ratio P for AR, HAR and RFSV predictors for $log(\sigma_{t+\Delta}^2)$</caption><tr><td>BCOM Δ = 1</td><td>0.846</td><td>0.838</td><td>0.805</td><td>0.83</td></tr><tr><td>BCOM Δ = 5</td><td>0.876</td><td>0.854</td><td>0.821</td><td>0.825</td></tr><tr><td>BCOM Δ = 21</td><td>0.956</td><td>0.937</td><td>0.874</td><td>0.854</td></tr><tr><td>INDU Δ = 1</td><td>0.451</td><td>0.446</td><td>0.444</td><td>0.458</td></tr><tr><td>INDU Δ = 5</td><td>0.617</td><td>0.612</td><td>0.532</td><td>0.541</td></tr><tr><td>INDU Δ = 21</td><td>0.858</td><td>0.857</td><td>0.716</td><td>0.699</td></tr><tr><td>USDEUR Δ = 1</td><td>0.530</td><td>0.514</td><td>0.507</td><td>0.521</td></tr><tr><td>USDEUR Δ = 5</td><td>0.611</td><td>0.581</td><td>0.532</td><td>0.544</td></tr><tr><td>USDEUR Δ = 21</td><td>0.755</td><td>0.728</td><td>0.618</td><td>0.638</td></tr><tr><td>IBOV Δ = 1</td><td>0.602</td><td>0.595</td><td>0.587</td><td>0.617</td></tr><tr><td>IBOV Δ = 5</td><td>0.779</td><td>0.754</td><td>0.68</td><td>0.691</td></tr><tr><td>IBOV Δ = 21</td><td>1.010</td><td>1.008</td><td>0.843</td><td>0.836</td></tr><tr><td>MICROSOFT Δ = 1</td><td>0.579</td><td>0.576</td><td>0.566</td><td>0.603</td></tr><tr><td>MICROSOFT Δ = 5</td><td>0.749</td><td>0.737</td><td>0.668</td><td>0.673</td></tr><tr><td>MICROSOFT Δ = 21</td><td>0.936</td><td>0.931</td><td>0.807</td><td>0.79</td></tr><tr><td>GOOGLE Δ = 1</td><td>0.500</td><td>0.497</td><td>0.492</td><td>0.529</td></tr><tr><td>GOOGLE Δ = 5</td><td>0.683</td><td>0.672</td><td>0.581</td><td>0.595</td></tr><tr><td>GOOGLE Δ = 21</td><td>0.864</td><td>0.861</td><td>0.729</td><td>0.722</td></tr><tr><td>SP400 Δ = 1</td><td>0.454</td><td>0.451</td><td>0.445</td><td>0.464</td></tr><tr><td>SP400 Δ = 5</td><td>0.616</td><td>0.601</td><td>0.525</td><td>0.538</td></tr><tr><td>SP400 Δ = 21</td><td>0.816</td><td>0.81</td><td>0.668</td><td>0.67</td></tr></table>

As we can see in Table 5.6, even though RFSV sometimes underperform AR, and HAR for $\Delta = 1$, it performs at least as good as the HAR when predicting more days ahead ($\Delta = 5, 21$) and outperforms the AR model.

Compared to AR, HAR whose parameters change through time, depend on the time horizon, need to be re-calibrated and even encounter calibration issues for some periods, the RFSV is more parsimonious since it only requires the parameter H to forecast the log-variance. In addition to that, the smoothness typically does not change over time or very slightly.

We notice that prediction through the RFSV can be linked to that of [67], where the issue of the prediction of the log-volatility in the multifractal random walk model of [12] is tackled. In this model,

$$ \mathbb{E}[\log(\sigma_{t+\Delta}^2) | \mathcal{F}_t] = \frac{1}{\pi} \sqrt{\Delta} \int_{-\infty}^{t} \frac{\log(\sigma_s^2)}{(t-s+\Delta)\sqrt{t-s}} dt, $$

which is the limit of our predictor when H tends to zero.

The prediction formula for the RFSV model can also be rewritten as

$$ \mathbb{E}[\log(\sigma_{t+\Delta}^2) | \mathcal{F}_t] = \frac{\cos(H\pi)}{\pi} \int_0^{+\infty} \frac{\log(\sigma_{t-\Delta u}^2)}{(u+1)u^{H+1/2}} du, $$

for a given small $\epsilon > 0$, let $r$ be the smallest real number such that

$$ \int_r^{+\infty} \frac{1}{(u+1)u^{H+1/2}} \le \epsilon. $$