I am struggling a bit with a coding project in R at the moment. I am to the price process of an arbitrary stock generated by a geometric brownian motion, then find the price for a European call option with arbitrary parameters and create a replicating portfolio for it. This first part has gone well.
The second part is essentially to repeat the process, but simulating the possibility of cracks in the market, which I have tried to do with a compound Poisson process, and generally it seems to have gone somewhat well. The plots for this process seem generally really close to what I want to simulate except for a couple of points.
Firstly, the interval on which the second simulation was done is different from the first one. For the second simulation the interval is (0, N=1460) increasing by increments of 1, whereas on the first one the interval is (0, T=1) increasing by increments of dt=T/N. Secondly the last period in my second process seems to be acting really erratic and producing incredibly strange values.
The second issue can at worst be fixed by ignoring the last period as the simulation of the replicating portfolio should already be close enough at the second to last period, but without the simulations running over the same time interval it is really hard for me to compare the two situations.
Please find the attached R-code
# Constants needed for computation
x0=100; mu=0.1; sigma=0.25; N=1460; K=100; T=1; dt=T/N; rho=0.05
# Seed and normal distribution
set.seed(1001)
X = rnorm(N)
X = append(X, 0, 0)
# Brownian motion
Bt = cumsum(X)*sqrt(dt)
# Date vector
t = seq(0, T, dt)
# 95% CI
lower = -1.96*sqrt(t); upper = 1.96*sqrt(t)
# Plotting the brownian motion and CI
plot(t, Bt, type="l", ylim=c(min(lower), max(upper)))
lines(t, lower)
lines(t, upper)
# Generating price process of the stock
xt = x0*exp((mu-sigma^2/2)*t+sigma*Bt)
# Expected value of stock at time t
ext = x0*exp(mu*t)
# 95% CI for stock price process
lxt = x0*exp((mu-sigma^2/2)*t-1.96*mu*sqrt(t))
uxt = x0*exp((mu-sigma^2/2)*t+1.96*mu*sqrt(t))
# Plotting price process and CI
plot(t, xt, type="l", ylim=c(min(lxt),max(uxt)))
lines(t, ext)
lines(t, lxt)
lines(t, uxt)
# Deriving the BSM formula
eta = sigma^(-1)/sqrt(T-t) * (log(xt/K)+rho*(T-t))
d1 = eta+sigma^2/2*sqrt(T-t)
# Finding the delta, amount of stocks held, at all times t (theta1)
Delta = pnorm(d1[1:N])
Delta = append(Delta, if(xt[N+1]>=K) 1 else 0)
# Value of stock position
stockVal = xt*Delta
# Price of option at time zero
p = x0*pnorm(d1[1]) - exp(-rho*T)*K*pnorm(d1[1]-sigma*sqrt(T))
# rlp = RiskLess Position
rlp = function(u) {
if (u == 1) return (p - Delta[1]*xt[1])
else return (rlp(u-1)*exp(rho*dt) - (Delta[u]-Delta[u-1])*xt[u])
}
# Plot the stock (black) and riskless (red) positions
pos = matrix(c(stockVal[-1], lapply(1:N,rlp)), ncol=2)
matplot(t[-1], pos, type="l", xlab="t", ylab="Positions", col=c(1,2))
# Construct replicating portfolio
rp = function(u) return(Delta[u]*xt[u] + rlp(u))
# Compare with buyer's position on EU call option
rp(N+1)
max(xt[N+1]-K,0)
rp(N+1) - max(xt[N+1]-K,0)
# Generate compound Poisson process
Nt <- rpois(N, lambda * T) # Number of jumps at each time step
Y <- numeric(N) # Initialize Y_t
for (i in 1:N) {
if (Nt[i] > 0) {
Y[i] <- sum(runif(Nt[i], min = -0.3, max = 0)) # Sum of jumps
} else {
Y[i] <- 0
}
}
# Brownian motion
Bt <- sqrt(dt) * cumsum(rnorm(N))
# Compute price process
Z <- (mu + 0.5 * sigma^2) * t + sigma * Bt + cumsum(Y)
S <- S0 * exp(Z)
# Plot
plot(S, type = 'l', xlab = 'Time', ylab = 'Price', main = 'Simulated Price Process')
I am generally quite inexperienced with coding, so it’s very hard for me to tell where exactly the code has gone wrong, but it seems to me that in the generation of the compound Poisson process, something must have happened.
Harald Straume is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.