Optimization of covered calls under uncertainty

Abstract

We present a two-stage stochastic program with recourse to construct covered call portfolios. To maximize the expected utility of a covered call portfolio, the model selects equity positions and call option overwriting weights for varying strike prices and expiry dates. Since the model has linear constraints and risk-averse utility functions are concave, the optimization problem is convex. The model is tested using 67 U.S. large-cap equities and optimizing the quadratic, negative exponential, and power utility functions. The expected utility is modeled as the average utility of the portfolio in a number of scenarios. Scenarios are first generated randomly then moment matching is employed, allowing the model to produce high quality results with a relatively small number of scenarios. To improve solution times we use a progressive hedging decomposition.

Appendices

Appendix 1

Here we describe the production of scenario prices $S_{1j}^i$ and $S_{2j}^{im}$. We produce returns using geometric Brownian motion with GARCH volatilities, then we adjust the resulting log-returns so that the sample volatilities match the volatilities implied by the market prices of ATM options, and so that the sample average simple returns matches the historical simple return.

We first assume that daily log-returns are jointly distributed according to a multivariate normal distribution with a constant vector of drifts, a constant correlation matrix, and a vector of variances each following a GARCH(1, 1) process. For a single asset j:

$$\begin{aligned} r(j,t,i)=\left( \mu _j- \frac{1}{2} \;\sigma (j,t,i)^2\right) +\sigma (j,t,i)Z(j,t,i) \end{aligned}$$

where

• r(j, t, i) is the log-return of asset j on day t in scenario i.

• $\mu _j$ is the historical drift of the daily log-return of asset j, gross of dividends

• $\sigma (j,t,i)^2$ is the variance of the daily log-return of asset j on day t in scenario i, modeled as a GARCH(1,1) process.

• Z(j, t, i) are standard normal random variables which are jointly distributed across assets: corr(Z(j, t, i), Z(k, t, i))=corr$(r(j,t,i),r(k,t,i))=\rho _{j,k}$.

The drifts and correlations are estimated using historical data from November 2006 to November 2016 and are assumed to be constant. GARCH(1,1) processes are fitted on the daily log-return residuals over the same time period. From fitting the GARCH processes we have variance estimates for all assets at time $t=1$, denoted by $\sigma (j)$. We use these volatilities as the starting point of the GARCH processes in all scenarios, i.e. $\sigma (j,t=1,i)=\sigma (j)$ for all scenarios i and all assets j. For each scenario i we can simulate $Z(j,t=1,i)$ for all assets j from a multivariate normal distribution with means equal to zero, variances equal to one, and covariance matrix $\rho =[\rho _{j,k}]$. Then, for all assets j in each scenario i, the log-return $r(j, t=1, i)$ can be computed and the GARCH variances can be updated through the equation $\sigma (j, t+1, i)^2 = \alpha ^0_j +\alpha _j \sigma (j,t,i)^2 Z(j, t, i)^2 + \beta _j \sigma (j,t,i)^2$, where $\alpha ^0_j$, $\alpha _j$, and $\beta _j$ are the fitted GARCH(1,1) parameters for asset j. We can then proceed to $t=2$. In this fashion log-returns are simulated one day at a time for $t=1$ to $T_1$.

For an asset j and a scenario i, the log-return during the first stage is given by:

$$\begin{aligned} r_{1j}^i=\sum _{t=1}^{T_1} r(j,t,i) \end{aligned}$$

Consider the sample expectation and variance of $r_{1j}^i$ across the simulated scenarios:

$$\begin{aligned} \text {var}(r_{1j})&\equiv {\hat{\sigma }}_j^2 T_1 \\ {\mathbb {E}}(r_{1j})\equiv M_j&\equiv \left( {\hat{\mu }}_j - \frac{1}{2}{\hat{\sigma }}_j^2\right) T_1 \end{aligned}$$

Suppose that we would like to transform the values of $r_{1j}^i$ so that the sample daily variance is equal to the daily variance implied by the market price of the ATM 2-week-to-maturity call option on asset j, ${\tilde{\sigma }}_j^2$, and so that the sample average simple return matches that implied by the historical rate $\mu _j$. Furthermore we would like to preserve higher order moments and correlations present in the samples. We can do this by applying the transformation:

$$\begin{aligned} {\tilde{r}}_{1j}^i&=\tau (r_{1j}^i-M_j)+{\tilde{M}}_j \\ \tau&= \frac{{\tilde{\sigma }}_j}{{\hat{\sigma }}_j}, \quad {\tilde{M}}_j=\left( {\mu }_j-\frac{1}{2}{\tilde{\sigma }}_j^2\right) T_1 \end{aligned}$$

Under this transformation it is straightforward to show:

$$\begin{aligned} {\mathbb {E}}(\text {exp}({\tilde{r}}_{1j}))&\approx \text {exp}({\mu }_j T_1) \\ \text {var}({\tilde{r}}_{1j})&={\tilde{\sigma }}_j^2 T_1 \\ \text {corr}({\tilde{r}}_{1j},{\tilde{r}}_{1k})&=\text {corr}({r}_{1j},{r}_{1k}) \\ {\bar{\mu }}_v({\tilde{r}}_{1j})&= {\bar{\mu }}_v({r}_{1j}) \quad v\in {\mathbb {N}} \end{aligned}$$

where ${\bar{\mu }}_v$ is the standardized moment of degree v.

We compute the price of asset j in scenario i at the end of the first stage: $S_{1j}^i=S_{0j}\text {exp}({\tilde{r}}_{1j}^i)-D_{1j}$, where if an ex-dividend date falls between $t=1$ to $T_1$ then $D_{1j}$ is the value of that dividend, or zero otherwise.

Production of scenarios for the second stage is nearly identical. In the second stage the GARCH volatilities are initialized to their final realized value from the first stage simulations, i.e. $\sigma (j,t=T_1+1,im)=\sigma (j,t=T_1+1,i)$ for all assets j and scenarios im. For the second stage the target variance ${\tilde{\sigma }}_j$ is set to be the inferred ATM implied volatility from $t=T_1$ to $t=T_1+T_2$ based on the daily volatilities implied by the ATM 4-week-to-maturity and 2-week-to-maturity options at time zero. This ensures that the volatility of log-returns from $t=0$ to $t=T_1+T_2$ matches the implied volatility of ATM 4-week-to-maturity options. The remainder of the procedure is unchanged, we set: $S_{2j}^{im}=S_{1j}^i\text {exp}({\tilde{r}}_{2j}^{im})-D_{2j}$.

Appendix 2

Here we describe the procedure used to estimate the market prices of call options at the beginning of the second stage in each scenario i, i.e. $C_{jl}^i$. Supposing that we have produced scenario prices $S_{1j}^i$ for all assets at the end of the first stage, we then require some method for estimating the market prices of call options such that they are consistent with the scenario asset prices. The maturity date is in $T_2$ days, i.e. at the end of the second stage. The risk-free rate is assumed to be known, and the assets’ spot prices are given by the scenario prices. From an option pricing perspective the only other parameters needed are the strike prices and volatilities.

We pick strike prices for each asset’s options by first assuming that the number of OTM options available for some asset is equal to the number of $T_1$-days-to-maturity OTM options which were available at the beginning of the first stage, i.e. $n_j^i=n_j^0$ for all assets j in each scenario i. Strike prices for a particular asset are generally given by multiples of some number centered around the spot price, e.g. multiples of 1, 5, 10, etc. We observe this strike price multiplier for each asset using the $n_j^0$ options available at time zero. We set the strike price of the ATM option on asset j at the beginning of the second stage in scenario i by rounding the spot price $S_{1j}^i$ down to the closest strike price multiple. Subsequent OTM strikes are set by increasing the ATM strike by the observed multiplier until $n_j^i=n_j^0$.

We require an estimate of the implied volatility curve where the expiry is in $T_2$ days. A substantial amount of literature exists in modelling implied volatility surfaces, much of which is collected in Gatheral (2006). In our case we may view the entire volatility surface for each asset at time zero. Since there are only $T_1$ days to the end of the first stage it is reasonable to assume that changes in the volatility surface over this time period ought to be small. We assume that for an asset j the implied volatility curve for $T_2$-days-to-maturity options at the beginning of the second stage, as a function of moneyness, is unchanged from that of time zero. This assumption is known as the ‘sticky delta’ rule. We estimate market prices by first computing implied volatilities of $T_2$-days-to-maturity options at time zero using a Cox-Ross-Rubinstein (CRR) binomial tree (see Cox et al. 1979 for details), then estimate options’ market prices at the beginning of the second stage using a CRR tree with the implied volatilities and the relevant strike prices. It is possible to apply sophisticated models if needed; since the scenario prices are exogenous inputs to the optimization model, prices generated using any model can easily be substituted.

Appendix 3

Here we prove that it is equivalent to optimize the expectation of two increasing utility functions which have the same Arrow–Pratt ARA coefficient.

Theorem

If two monotonically increasing functions $f:{\mathbb {R}}\mapsto {\mathbb {R}}$ and $g:{\mathbb {R}}\mapsto {\mathbb {R}}$ satisfy the condition $f''(x)/f'(x) = g''(x)/g'(x)$ then $\exists \; a,b \in {\mathbb {R}}$ such that $f(x)=a+b\,g(x)$.

Proof

Note that $f''(x)/f'(x) = \frac{\text {d}}{\text {d}x} \text {log}(f'(x))$. Then:

$$\begin{aligned} \frac{f''(x)}{f'(x)}&=\frac{g''(x)}{g'(x)} \\ \frac{\text {d}}{\text {d}x} \text {log}(f'(x))&= \frac{\text {d}}{\text {d}x} \text {log}(g'(x)) \\ \text {log}(f'(x))&= \text {log}(g'(x)) + c \\ f'(x)&= \text {e}^c \, g'(x) \\ f(x)&= a+b\, g(x) \end{aligned}$$

where $c \in {\mathbb {R}}$ is a constant and defining $b=\text {e}^c$. $\square$

Applying the theorem above, two monotonically increasing utility functions with the same ARA coefficients differ by only a linear transformation. Maximizing ${\mathbb {E}}(U(x))$ is equivalent to maximizing ${\mathbb {E}}(a+b\, U(x))$ given that $b>0$ since both functions are increasing.

Appendix 4

Here we demonstrate that setting per element values of r in the progressive hedging algorithm is equivalent to using a single value of r for all variables and rescaling the decision variables.

Theorem

When using the progressive hedging algorithm with a penalty parameter r rescaling a decision variable $x_k$ by multiplying it by a constant $n_k$ is equivalent to using the variable specific penalty value $r_k =r\, n_k^2$.

Proof

Consider the penalty associated with a single decision element $x_k$ in iteration $\nu$ of the progressive hedging algorithm for some particular subproblem:

$$\begin{aligned}&V_{\nu , k} x_k + 0.5 r (x_k - {\hat{x}}_{\nu , k})^2 \\&\quad = \left( \sum _2^\nu r(x_{\nu , k} - {\hat{x}}_{\nu , k}) \right) x_k + 0.5 r (x_k - {\hat{x}}_{\nu , k})^2 \end{aligned}$$

Consider replacing $x_k$ with a rescaled variable $y_k=n_k x_k$. The optimization problem is unchanged assuming we make the substition $x_k = y_k / n_k$ in the subproblem’s constraints and original objective function. Consider the penalty that progressive hedging would impose on $y_k$ in iteration $\nu$:

$$\begin{aligned}&\left( \sum _2^\nu r(y_{\nu , k} - {\hat{y}}_{\nu , k}) \right) y_k + 0.5 r (y_k - {\hat{y}}_{\nu , k})^2 \end{aligned}$$

Transforming back to the original variable (assuming for now that ${y}_{\nu ,k}=n_k {x}_{\nu , k}$, which also implies ${\hat{y}}_{\nu ,k}=n_k {\hat{x}}_{\nu , k}$):

$$\begin{aligned} =&\left( \sum _2^\nu r(n_k x_{\nu , k} - n_k {\hat{x}}_{\nu , k}) \right) n_k x_k + 0.5 r (n_k x_k - n_k {\hat{x}}_{\nu , k})^2 \\ =&\left( \sum _2^\nu r\, n_k^2(x_{\nu , k} - {\hat{x}}_{\nu , k}) \right) x_k + 0.5 r\,n_k^2 (x_k - {\hat{x}}_{\nu , k})^2 \\ =&\left( \sum _2^\nu r_k(x_{\nu , k} - {\hat{x}}_{\nu , k}) \right) x_k + 0.5 r_k (x_k - {\hat{x}}_{\nu , k})^2 \end{aligned}$$

where we defined $r_k=r\, n_k^2$. This is identical to the original penalty for $x_k$ except with r replaced by $r_k$. It remains to show ${y}_{\nu ,k}=n_k {x}_{\nu , k}$. The first iteration is solved without penalties, thus for each subproblem ${y}_{2,k}=n_k {x}_{2, k}$ is true, and ${\hat{y}}_{2,k}=n_k {\hat{x}}_{2, k}$ is true. Suppose that we use the penalty parameter $r_k=r\, n_k^2$ for element $x_k$ in the original subproblems and the penalty r in the rescaled subproblems. Then the penalties are identical in the second iteration since ${y}_{2,k}=n_k {x}_{2, k}$, and the problems have the same solution, ${y}_{3,k}=n_k {x}_{3, k}$. By induction this is true for all subsequent iterations. Therefore, the progressive hedging algorithm with $x_k$ rescaled by $n_k$ is equivalent to performing progressive hedging using the value $r_k=r\, n_k^2$ for element $x_k$. $\square$

Since the progressive hedging penalty is separable across variables, the proof holds when multiple variables are rescaled.