Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: THEend8_
CS 335 Assignment
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
∆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
Si+1 = Si + Si(r∆t+ σφi
√
∆t), (1)
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
V (S0, 0) = e
−rT
(∑M
m=1 Payoff((SN )m)
M
)
.
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 0.08
Time to expiry T 1.0 years
Strike Price K $101
Initial asset price S0 $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
Π(0) = −V (S(0), 0) + α(0)S(0) +B(0) = 0
α(0) =
∂V
∂S
(S(0), 0)
B(0) = V (S(0), 0)− α(0)S(0).
The hedge is rebalanced at discrete times ti. Defining
αi =
∂V
∂S
(Si, ti)
Vi = 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
Πi = −Vi + αiSi +Bi.
2
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(ti+1) = S(ti) exp
[
(µ− σ2/2)(ti+1 − ti) + σφ
√
ti+1 − ti
]
, (2)
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 No
Volatility σ 0.16
Risk free rate r 0.06
Time to expiry T 2 years
Strike Price K $100
Initial asset price S0 $100
Real world drift rate µ 0.08
Payoff type Put
Number of Hedging rebalances 10, 100, 1000, 10000
Table 2: Data for Hedging Simulations (Q2).
3. (12 marks) CPPI simulation.
3
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 1/250
initial cash 100
initial risky asset position α0 0
initial risky asset price S0 100
Number of simulations 80000
Table 3: Data used in the CPPI strategy experiments.
Let
B(ti) = Bi = Amount in risk free asset at ti
S(ti) = Si = Price of risky asset at ti
α(ti) = αi = Number of units of risky asset at ti
Π(ti) = Πi = Bi + αiSi = Total portfolio at ti
M(ti) = Mi = CPPI multiplier at ti
F (ti) = Fi = 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
[
max(0, Bie
r∆t + αiSi+1 − Fi+1)
Si+1
]
(3)
Bi+1 = Bie
r∆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 No
Volatility σ 0.35
Risk free rate r 0.07
Time to expiry T 1.5
Strike price K 101, 115, 130
Initial asset price S0 100
Initial investment Sinit 100
Real world drift µ 0.09
Payoff Call
Number of simulations 80000
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.