Optimizing a Multi-Asset Portfolio of Long Call Options

# Call options portfolio: what would the underlyings' moments to be maximized?


# Call options portfolio: what would the underlyings' moments to be maximized?












Let you have only three underlyings, like SPY, TLT and GLD, and you want to buy $n_{1}$ Call options on SPY, $n_{2}$ Call options on TLT and $n_{3}$ Call options on GLD... with a limited budget, that is obvious.

You know those three underlyings, $S_{i}$ where $i=1,2,3$, are correlated as explained by $P$ correlation matrix and $V$ annualized covariance matrix (in case, with implied volatilities in main diagonal).

Being moneyness and maturity equal for each Call option written on each underlying, my goal is to find the value of $n_{1}, n_{2}$ and $n_{3}$ in order to make often good profits over time.

What was my step-by-step proceedings in such a case?

- Calculating pseudo-square root matrix of $V$, that is $\Sigma$, via eigenvalue decomposition;

- setting an horizon of 30 days, that is $T=30/365=0.082$ years;

- simulating 10,000 correlated multi-dimensional GBM stochastic processes with risk free rate drift, $\Sigma$ diffusion matrix and $T$ terminal time;

- calculating arithmetic returns separately for each underlying as $R_{T}=\frac{S_{T}}{S_{0}}-1$.

If I had bought the underlying securities, I would have linear P&L and I could choose my favorite performance (risk) measure to be maximized (minimized), like Omega, Expected Shortfall, Tail Dependance, Drawdown etc.; but in this case I'm dealing with options, which have not linear P&L due to convexity (Gamma).

This means you're looking for realized volatility (dispersion of $S_{T}$) higher than implied volatility regardless of weighted returns trend, because Call options already defend you against drawdowns and tail risk.

Hence my question is: making $n_{1}, n_{2}$ and $n_{3}$ varying over time, what's the measure I should maximize or minimize to

- maximize my returns;

- maximize my Sharpe;

- making my return shape distribution the best possible

...buying only Call options on those underlying securities?

Here is the `R` reproducible code to maximize long Call options modified Sharpe ratio. In the real world book's options prices should be used, but in this example I've used BMS formula with GARCH(1,1) volatility to estimate them and their change over time.

Please, tell me if you think there are theoretical and/or practical mistakes:

```
# ======================================================== #
# Managing options portfolio risk                          #
# Stochastic optimization of several correlated securities #
# By Lisa Ann for Quantitative Finance @ StackExchange                                                  #
# ======================================================== #

# Contents:

# 1. Loading packages
# 2. Downloading datas and preparing environment
# 3. Estimating orthogonal GARCH(1,1) covariance matrix
# 4. Simulating stochastic n-dimensional diffusion processes with 'n' equations
# 5. Correlated geometric Brownian motions returns
# 6. Pricing options at the end of simulation
# 7. Optimization of portfolio modified Sharpe

# *********************************
# 1. Loading packages
# *********************************

library(DEoptim)
library(quantmod)
library(RQuantLib)
library(Rsolnp)
library(rugarch)

# 'yuima' package from R-Forge

install.packages("yuima", repos="http://R-Forge.R-project.org")
library(yuima)

# You can download 'rmgarch' package from:
# https://rmgarch.googlecode.com/files/rmgarch_1.00-4.zip

library(rmgarch)

# *********************************
# 2. Downloading datas and preparing environment
# *********************************

# We will start with few assets which represent main asset classes: SPY for 
# equity, TLT for Government bonds and GLD for gold. We will add more assets
# later on

tickers <- c('SPY', 'TLT', 'GLD')
getSymbols(tickers, from = '1950-01-01')

# 3-Month London Interbank Offered Rate (LIBOR), based on U.S. Dollar

getSymbols('USD3MTD156N', src = 'FRED')

# Merging close-to-close arithmetic returns

xinit <- as.vector(tail(na.omit(merge(Cl(SPY), Cl(TLT), Cl(GLD))), 1))
data <- na.omit(merge(ClCl(SPY), ClCl(TLT), ClCl(GLD)))
xyplot(data)

# Risk free rate proxy

r.f <- as.numeric(tail(USD3MTD156N, 1) / 100)

# *********************************
# 3. Estimating orthogonal GARCH(1,1) covariance matrix
# *********************************

# We have to set our time horizon

n.ahead <- 30

# Covariance matrix forecast

spec <- gogarchspec(mean.model = list(model = 'constant'), 
                    variance.model = list(model = 'apARCH',
                                          garchOrder = c(1,1)),
                    distribution.model = 'manig')
fit <- gogarchfit(spec = spec, data = data)
forecast <- gogarchforecast(fit = fit, n.ahead = n.ahead)
V <- rcov(forecast)[[1]][,,n.ahead]
colnames(V) <- rownames(V) <- tickers
V

# *********************************
# 4. Simulating stochastic n-dimensional diffusion processes with 'n' equations
# *********************************

# Eigenvalue decomposition of coavariance matrix:
# a is the matrix for which we want to square root, then
# a.eig <- eigen(a)
# a.sqrt <- a.eig$vectors %*% diag(sqrt(a.eig$values)) %*% solve(a.eig$vectors)

s.eig <- eigen(V)
sigma.ch <- s.eig$vectors %*% diag(sqrt(s.eig$values)) %*% t(s.eig$vectors)
colnames(sigma.ch) <- rownames(sigma.ch) <- tickers
sigma.ch

# Stochastic differential equations

n <- ncol(data)
drift <- paste('r.f * x', 1:n, sep = '')
diffusion <- matrix(NA, n, n)

for(i in 1:n){
  for(j in 1:n){
    diffusion[i,j] <- paste(sigma.ch[i,j], ' * x', j, sep = '')
  }
}

object <- setModel(drift = drift, diffusion = diffusion, 
         solve.variable = paste('x', 1:n, sep = ''))

# Sampling grid properties: 'n.ahead' days and number of intervals

sampling <- setSampling(Terminal = n.ahead / 365, n = n.ahead * 1)

# Plotting simulations to check SDEs behaviour follows covariance matrix

plot(simulate(object = object, xinit = xinit, sampling = sampling), lwd = 2,
     type = 'o', col = 1:n)

# *********************************
# 5. Correlated geometric Brownian motions
# *********************************

# Starting values

S.0 <- matrix(xinit, nrow = 10000, ncol = n, byrow = TRUE)

# Terminal values

S.T <- matrix(NA, nrow = 10000, ncol = n)

for(i in 1:10000){
  for(j in 1:n){
    sim <- simulate(object = object, xinit = xinit, sampling = sampling)
    S.T[i,j] <- tail(sim@data@zoo.data[[j]], 1)
  }
}

colnames(S.T) <- tickers
par(mfrow = c(n,n))

for(j in 1:n){
  hist(S.T[,j], xlab = tickers[j], freq = FALSE,
       main = paste('Probability of', tickers[j], 'being at', 
                    expression(S(T)), n.ahead, 'periods from today'))
  lines(density(S.T[,j]), lwd = 2)
}

par(mfrow = c(1,1))

# *********************************
# 6. Pricing options at the end of simulation
# *********************************

# Now we have to price Call options at time 'T' to see true performances
# Estimating volatility of underlyings

spec <- ugarchspec(variance.model = list(model = 'sGARCH', garchOrder = c(1,1)),
                   mean.model = list(armaOrder = c(0,0)),
                   distribution.model = 'norm')

# Function to estimate underlying volatility

hv <- function(x, t, n.ahead){
  fit <- ugarchfit(spec = spec, data = x)
  forc <- ugarchforecast(fit, n.ahead = n.ahead)
  sigma <- forc@forecast$sigmaFor[t]
  return(sigma)
}

# Volatilities at '0' and at 'T'

sigma.0 <- apply(X = data, MARGIN = 2, FUN = hv, t = 1, 
                 n.ahead = n.ahead) * sqrt(252)
sigma.0.vec <- matrix(sigma.0, nrow = nrow(S.0), ncol = n, byrow = TRUE)
sigma.T <- apply(X = data, MARGIN = 2, FUN = hv, t = n.ahead, 
                 n.ahead = n.ahead) * sqrt(252)
sigma.T.vec <- matrix(sigma.T, nrow = nrow(S.T), ncol = n, byrow = TRUE)
barplot(rbind(sigma.0, sigma.T), beside = TRUE, legend.text = TRUE,
        ylab = 'Annualized Volatility')

# Options prices at 'S.0'
# We consider (ATM + 1) strikes for Call options

strike <- ceiling(S.0) + 1
call.0 <- matrix(NA, nrow = nrow(S.0), ncol = n)

for(i in 1:nrow(S.0)){
  for(j in 1:n){
    call.0[i,j] <- EuropeanOption(type = 'call', underlying = S.0[i,j], 
                                  strike = strike[i,j], dividendYield = r.f, 
                                  riskFreeRate = r.f, 
                                  maturity = (n.ahead * 2) / 365,
                                  volatility = sigma.0.vec[i,j])$value  
  }
}

# Options prices at 'S.T'

call.T <- matrix(NA, nrow = nrow(S.T), ncol = n)

for(i in 1:nrow(S.T)){
  for(j in 1:n){
    call.T[i,j] <- EuropeanOption(type = 'call', underlying = S.T[i,j], 
                                  strike = strike[i,j], dividendYield = r.f, 
                                  riskFreeRate = r.f, maturity = n.ahead / 365,
                                  volatility = sigma.T.vec[i,j])$value    
  }
}

# Call buyer returns

call.r <- call.T / call.0 - 1

# *********************************
# 7. Optimization of portfolio modified Sharpe
# *********************************

# Call weighted returns function

call.w <- function(w){
  ret <- call.r %*% w
  mu <- mean(ret)
  q <- quantile(ret, probs = .05)
  risk <- mean(ret[ret < q])
  return(- (mu - r.f) / risk)
}

# Optimization

par <- DEoptim(fn = call.w, lower = rep(10, n), upper = rep(100, n),
               control = DEoptim.control(strategy = 6, 
                                         itermax = 1000))$optim$bestmem
par / sum(par)
round(solnp(pars = runif(3, 0, 1), fun = call.w, eqfun = function(w){sum(w)},
            eqB = 1, LB = rep(0, n), UB = rep(1, n), 
            control = list(innter.iter = 10000))$pars, 4)
```

At the end of optimization, my budget is allocated on SPY, TLT and GLD according to optmal weights.