CS 335 - Winter 2023: Assignment 4
Submit all components of your solutions (written/analytical work, code/scripts/notebooks, figures, plots, etc.) to CrowdMark in PDF form in the section for each question.
You must also separately submit a single zip file containing any and all code/notebooks/scripts you write to the DropBox on LEARN in runnable format (that is, the files in the zip should be of type .ipynb).
For full marks, you must show your work and adequately explain your logic!
1. (10 marks) Monte Carlo option pricing.
As we derived in class, for the purposes of pricing options, we can pretend that the asset price
S evolves in the risk-neutral world as
dS = rSdt + σ(S, t)SdZ,
where r is the risk free rate, σ(S,t) is the volatility (as a function of S and t), and dZ is the increment of a Weiner process.
Let the expiry time of an option be T, and let N=T,
Si = S(i?t),
where ?t is the time step size. Then, given an initial price S0, M different simulations of the path of the asset should be generated using the simple time-stepping rule
where φi ～ N(0,1) at each step. If the value of SN = S(T) generated by the mth simulation is denoted by (SN)m, then an approximate value (i.e., Monte Carlo estimate) of the option price is given by
rT PMm=1 Payoff((SN)m)! M .
V (S0, 0) = e
For the volatility function, we will use σ(S, t) = γ/S where γ = 20 is a constant.
Since the basic time stepping method of Equation (1) could potentially produce negative values, for this question we will simply set the asset price to ε = 10?8 whenever timestepping produces Si+1 ≤ 0.
Prepare a Jupyter notebook that implements this method in Python/NumPy. Use the data for computing European puts given in Table 1. For a fixed ?t, write code to plot the computed option value versus number of simulations M, for each of
M = 1000, 2000, 4000, 8000, 16000, 32000, 64000, 96000, 128000, 192000, 256000. 1
Generate such a plot for the following values of ?t = 10/250, 2.5/250, 1/250. (We often assume 250 days in a financial year). Title each of your 3 plots with the corresponding timestep size.
Comment on the behavior you observe from these plots.
Additionally, for the smallest ?t value, write code to generate a simple text table where each row contains M, the estimated option value, and the lower and upper bound values of the 95% confidence interval, for that value of M.
Be sure to vectorize the inner loop of your code. That is, for a given choice of ?t and M, there should only be one explicit loop (i.e., over the time steps) so that all of the M simulations are computed simultaneously with vector operations, rather than one entry at a time.
Risk free rate r Time to expiry T Strike Price K Initial asset price S0
0.08 1.0 years $101 $100
γ 20 ε 10?8
Table 1: Parameters for Q1.
2. (10 marks) Delta hedging: single stochastic paths.
As we have seen in lecture, discrete delta hedging consists of a portfolio Π(t) with components:
A short option position ?V (S(t), t).
Long α(t) shares at price S(t).
An amount in a risk-free bank account B(t).
Initially, we have
The hedge is rebalanced at discrete times ti. Defining
αi = ?V(Si,ti)
then we have to update the hedge by purchasing αi ? αi?1 shares at t = ti, so that updating our share position requires
S(ti) (αi ? αi?1)
in cash, which we borrow from the bank. At the instant after the rebalancing time ti, the value of the portfolio is
Develop Python/NumPy code in a Jupyter notebook for simulating a stochastic path of a delta hedge. At the end of each discrete hedging interval, update your portfolio, using the current asset price, and the corresponding values of option and delta (using blsprice, blsdelta). (Note that you do not need to rebalance at the final time, t = T .)
Compute the relative hedging error at expiry time t = T , i.e. Relative Hedge Error = e?rT Π(T ) .
|V (0)|
Your code should be able to produce plots showing: the stock price Si, the value in the risk free account (Bi), the stock holding (αi Si), and the total portfolio value Πi at each point along this single stochastic path. The plots should be loosely similar to Figure 8.1.3 in the course notes. (Of course, your actual values will be different, due to the randomness inherent in the algorithm). Use the PyPlot function legend to identify the different curves.
Using the code above, generate four different plots, using 10, 100, 1000, 10000 timesteps / inter- vals in [0,T].
Generate an additional plot of number of rebalances vs. the corresponding absolute value of the relative hedge error. That is, it should have rebalance counts on the x-axis and error on the y-axis.
To evolve the stock price S for this question, we will assume GBM and use the following exact update rule for GBM:
S(t ) = S(t ) exp h(μ ? σ2/2)(t ? t ) + σφpt ? t i , (2) i+1 i i+1 i i+1 i
where φ ～ N(0,1).
Use the data in Table 2. There is no need to vectorize this code. Comment on what you observe.
Early Exercise Volatility σ
Risk free rate r
Time to expiry T
Strike Price K
Initial asset price S0
Real world drift rate μ Payoff type
Number of Hedging rebalances
T (investment horizon) 2.0 σ (volatility) 0.32 μ (real world drift) 0.12
P0 (initial wealth) 100 r (risk free rate) 0.06
Rebalancing interval
initial cash 100
initial risky asset position α0 0 initial risky asset price S0 100
Number of simulations
Data used in the CPPI strategy experiments.
Amount in risk free asset at ti Price of risky asset at ti
Number of units of risky asset at ti
Bi + αiSi = Total portfolio at ti CPPI multiplier at ti
CPPI floor at ti
A CPPI strategy is determined by specifying a floor Fi and a multiplier Mi. The state of the strategy is given by αi and Bi. At each rebalancing date ti+1, the following reallocation takes place
αi+1 = Mi+1 S (3)
i+1 Bi+1 = Bier?t ? (αi+1 ? αi) Si+1.
Prepare a Jupyter notebook that simulates CPPI, assuming GBM, using Python/NumPy. Use the exact time stepping method of Equation (2) for GBM, and the data in Table 3.
Measure investment performance in terms of R = log(Π(T)/Π(0)). Run a simulation for each of the following parameter pairs: (F, M ) = (0, 1), (0, 0.5), (0, 2), (85, 2), (85, 4), and generate a table containing F , M , and the mean, standard deviation, 95% VAR, and 95% cVAR of R. That is, F and M are considered to be constants in this case (rather than time-varying functions). Note also that the VAR and cVAR values will necessarily be approximations based on the data you generate - to estimate them, consider first sorting the vector of R results you compute.
For each of the above cases, use matplotlib’s hist function to generate a plot of the (approx- imate) probability density function of R, using 200 bins and setting density=True. Title each plot with the corresponding F and M values, e.g., “F = 0, M = 1”.
Note that initially you start out with cash and no stock, and the first rebalance occurs at t = 0+, i.e., determine α0 and B0 from equation (3) with α?1 = 0, B?1 = initial cash, ?t = 0.
Be sure to vectorize your code. Comment on what you observe from the plots and table. 4
Early exercise Volatility σ
Risk free rate r Time to expiry T Strike price K Initial asset price S0 Initial investment Sinit Real world drift μ Payoff Number of simulations
Table 4: Data for Covered Call Simulations (Q4).
4. (8 marks) Covered call writing.
A covered call strategy works as follows:
Assume we have an initial investment of Sinit, which is equal to the price of one share.
At t = 0, buy one share using Sinit, and sell one call option with expiry time T and strike K (i.e., we are now short the call). Put the cash received from selling the call into a risk-free bank account to earn interest.
At t = T, sell the stock, pay off the short position in the call, and collect the amount in the risk free account, which gives your final value in cash.
Assume that the price (at t = 0) of the call option is determined by the Black-Scholes model (i.e., using blsprice).
Write Python/Numpy code in a Jupyter notebook to simulate this strategy for three values of the strike K as given in Table 4. Assume the asset follows GBM (under the real world measure), i.e.,
dS = μSdt + σSdZ (4)
where as usual σ is the volatility, μ is the real world drift rate, and dZ is the increment of a Weiner process. Since σ and μ are assumed constant, we can again use the exact GBM solution in our simulations.
Measure investment performance using the value
log(final value/Sinit)
i.e., the log return on initial investment. Generate a table showing mean, standard deviation, 95% VAR, and 95% cVAR corresponding to each of the three strike values K of the call as listed in the table.
For the case with the largest strike price, generate a plot of the (approximate) probability density function of the log return. Use 80 bins.
Comment on what you observe.