%% Analyzing Investment Strategies with CVaR Portfolio Optimization in MATLAB - Generate Scenarios
%
% Robert Taylor
% The MathWorks, Inc.
% Copyright (C) 2012 The MathWorks, Inc.
%% Introduction
% This script generates scenarios for a covered-call strategy with a portfolio of stocks. The stocks
% are assumed to have no dividends (dividends have been incorporated into total return prices in the
% data) and the call options are assumed to be American options. The covered-call strategy has
% scenarios for uncovered stocks and for covered positions at two different strike prices.
% Subsequent optimization steps can use any combination of these three possible collections of
% scenarios. The scenarios are written to the file |BuyWriteScenarios.mat|.
%% Get Stock Total Return Price Data
% The example starts with data from the file |BuyWriteStockUniverse.mat| that contains a fints
% object, Universe, with asset symbols and 15 years of daily total return prices for 26 stocks.
load BuyWriteStockUniverse Universe
% Extract information from fints object Universe.
Assets = fieldnames(Universe, 1); % stock symbols
X = fts2mat(Universe); % stock total return prices
numassets = numel(Assets); % number of stocks in universe
%% Calibrate Geometric Brownian Motion
% Under an assumption that stock prices are geometric Brownian motion processes, the next step
% calibrates a continuous multivariate geometric Brownian motion process based on discrete samples
% of total return prices.
t0 = busdate(Universe.dates(1), -1); % first business day prior to first date in data
T = Universe.dates(end); % last date in data
% Determine average number of days in a year for the historical period. In this example, the period
% is known and fixed to be 15 years.
average_days_per_year = (T - t0)/15;
% Form "dates" for data relative to t0 in years based on average number of days in a year.
t = linspace(0, (T - t0), 1 + size(X,1))'/average_days_per_year; % "dates" from t0 for data
% Calibrate the model with maximum likelihood estimation.
tic
[x0, mu, gamma] = gbm_calibration(t(1), X, t(2:end));
toc_cal = toc;
fprintf('Time to calibrate parameters %g\n',toc_cal);
% Form various ancillary statistics from the estimated diffusion matrix.
covar = gamma*gamma';
[sigma, rho] = cov2corr(covar);
sigma = sigma(:);
% List calibrated parameters.
fprintf('Calibrated multivariate geometric Brownian motion parameters ...\n');
fprintf(' %5s %10s %10s %10s\n','Asset','Initial','Drift','Diffusion');
fprintf(' %5s %10s %10s %10s\n','','Price','Parameter','Parameter');
for i = 1:numassets
fprintf(' %5s %10.2f %10.2f %10.2f\n',Assets{i},x0(i),100*mu(i),100*sigma(i));
end
%% Set up Investment Period for Scenarios
% To generate scenarios, it is necessary to specify a fixed period of time called an investment
% period. During an investment period, numerous actions can take place that are independent of the
% specific features of a particular financial instrument. To set up an investment period, pick start
% and end dates. For this strategy a good choice is to select a single month.
startdate = datenum('28-Sep-2012'); % start date for investment period
enddate = datenum('31-Oct-2012'); % end date for investment period
numdays = daysdif(startdate, enddate, 13); % number of trading days in an investment period
t0 = 0; % set initial time to be 0 years
T = numdays/252; % terminal time is T years
%% Simulate Geometric Brownian Motion over an Investment Period
% Given an investment period and a calibrated geometric Brownian motion process, simulate stock
% prices for the underlying assets. Since the ultimate goal is to create a collection of scenarios
% for portfolio optimization, the number of scenarios is the number of trials for the stochastic
% process (note that a single trial has a path for each asset at the periodicity of the simulation
% and that the total return for the path is the scenario).
% Although the sampling rate for the simulation can be set to any value, this example uses half-hour
% increments. For the US markets, a trading day has 13 half-hour increments so this is the number of
% "ticks" per day.
num_ticks_per_day = 13; % number of "tick" periods per day
% To simulate a multivariate geometric Brownian motion, three control parameters must be set that
% are the number of trials, the number of samples per trial, and the number of time steps between
% samples. Note that the number of trials is the number of scenarios to be generated for CVaR
% portfolio optimization and a good "typical" number is 20000.
numtrials = 20000; % number of scenarios/trials/paths/realizations
numsamples = numdays*num_ticks_per_day; % number of samples/observables per scenario
numsteps = 1; % number of steps between samples
dt = (T - t0)/numsamples; % time interval between samples in years
% Simulate multivariate stochastic process.
% Xsde is a gbm object with the parameters from the calibration step for the process.
% X is a (1 + numsamples) by numassets by numtrials array of sample paths for the process.
% t is a (1 + numsamples) vector that contains the "times" for each sample in the process.
tic
Xsde = gbm(diag(mu), diag(sigma), 'starttime', t0, 'startstate', x0, 'correlation', rho);
[X, t] = Xsde.simulate(numsamples, 'deltatime', dt, 'ntrials', numtrials, ...
'nsteps', numsteps, 'antithetic', true);
toc_gbm = toc;
fprintf('Time to simulate gbm process %g with %g samples\n', ...
toc_gbm,(1 + numsamples)*numassets*numtrials);
%% Set up Covered-Call Strategy
% To gather the necessary inputs to generate scenarios for a covered-call strategy, the final step
% requires specification of various parameters that control the simulation. The two parameters that
% are of primary focus are initiating, which specifies whether initial call premia are to be
% included in scenario returns, and exercise_likelihood, which specifies the characteristics of the
% likelihood of early exercise. Although some guidelines to set these parameters are provided here,
% the function |covered_engine.m| contains a complete description of each of these parameters.
% Fund features.
initial_equity = 1.0e6; % initial wealth invested in portfolio
distribution = 0; % annualized fund distribution rate
risk_free_rate = 0.0015; % Treausury bill rate (annual yield), proxy for risk-free rate
no_reinvestment = false; % true if no re-investment after assignment, false otherwise
initiating = true; % true if initiating a call position, false if an ongoing position
% Call option features.
% Use a fixed expiration date for calls, including any calls that are rewritten. If the calls expire
% during the investment period, shift to the next contract expiration date. Select a date after the
% start date of the investment period. For this example with investment period end date, possible
% call expiration dates are: 19-Oct-2012, 16-Nov-2012, 21-Dec-2012, 18-Jan-2013, 15-Feb-2013,
% 15-Mar-2013, ... .
% The exercise_likelihood parameter is either NaN (to use a heuristic model) or a probability of
% early exercise if stock price >= strike price between start of investment period and option
% expiration.
near_contract = datenum('21-Dec-2012'); % near contract call option expiration date
next_contract = datenum('15-Mar-2013'); % next contract call option expiration date
contract_expiration = daysdif(startdate, near_contract, 13)/252; % "anchor" expiration dates
next_contract_expiration = daysdif(startdate, next_contract, 13)/252;
strike_cushion = 0.05; % "cushion" to set strike price above current stock price
exercise_likelihood = 0.99; % likelihood of early exercise if stock price >= strike price
% Costs/frictions.
% For the simulation, basic transaction costs are 5 c/share for stocks and 20 c/contract for
% options. Average opportunity costs are assumed to be 1 c/share for stocks and 2 c/contract for
% options.
stock_cost = 0.06; % 6 c/share
contract_cost = 0.22; % 22 c/option contract
% Assignment and re-investment delays.
confirmation_delay = 2; % 1 hour delay to confirm option assignment
reinvestment_delay = num_ticks_per_day; % 1 day to reinvest after confirmation of assignment
settlement_delay = 3*num_ticks_per_day; % 3 days to settle after reinvestment
%% Generate Scenarios for Uncovered and Covered Positions
% Generate scenarios by looping over the number of scenarios and the number of assets. The function
% |covered_engine.m| simulates covered-call scenarios and the function |uncovered_engine| simulates
% uncovered scenarios. Note that subsequent portfolio optimization steps may allocate weights to
% mixtures of covered and uncovered positions to form partial coverage of positions.
% WARNING: This step is slow.
tic
fprintf('Generating scenarios ...\n');
XU = zeros(numtrials, numassets);
XC1 = zeros(numtrials, numassets);
XC2 = zeros(numtrials, numassets);
for iter = 1:numtrials
for i = 1:numassets
rU = uncovered_engine(X(:,i,iter), T, ...
initial_equity, distribution, risk_free_rate, stock_cost);
rC1 = covered_engine(X(:,i,iter), T, mu(i), sigma(i), ...
initial_equity, distribution, no_reinvestment, initiating, ...
strike_cushion, contract_expiration, next_contract_expiration, risk_free_rate, ...
stock_cost, contract_cost, exercise_likelihood, ...
confirmation_delay, reinvestment_delay, settlement_delay);
rC2 = covered_engine(X(:,i,iter), T, mu(i), sigma(i), ...
initial_equity, distribution, no_reinvestment, initiating, ...
2*strike_cushion, contract_expiration, next_contract_expiration, risk_free_rate, ...
stock_cost, contract_cost, exercise_likelihood, ...
confirmation_delay, reinvestment_delay, settlement_delay);
XU(iter, i) = rU;
XC1(iter, i) = rC1;
XC2(iter, i) = rC2;
end
if mod(iter, 1000) == 0
fprintf('Completed %d scenarios ...\n',iter);
end
end
toc_scenarios = toc;
fprintf('Time to simulate scenarios %g\n',toc_scenarios);
%% Save scenarios
% Save scenarios in the file |BuyWriteScenarios.mat|.
save('./data/BuyWriteScenarios.mat','Assets','XU','XC1','XC2','exercise_likelihood','initiating');