Question 1

Here in question, we use Euler schemes to approximate the stochastic process

Before processing, noticing R has bugs in calculating nth root, define a function nthroot to avoid NA in calculating negative value’s nth root

nthroot = function(x,n) {
  (abs(x)^(1/n))*sign(x)
}
# Define SDE
sdeGeneratorX = function(N=10^4,T=1,x0,theta=c()){
  #set.seed(123)
  dt = T/N
  X = numeric(N+1)
  X[1] =x0
  W = rnorm(N)
  for (i in 1:N) {
    X[i+1] = X[i] + (theta[1]+theta[2]*X[i])*dt + theta[3]*sqrt(dt)*W[i]
  }
  return(X[N+1])
}
sdeGeneratorY = function(N=10^4,T=1,y0){
  #set.seed(123)
  dt = T/N
  Y = numeric(N+1)
  Y[1] =y0
  W = rnorm(N)
    for (i in 1:N) {
      Y[i+1] = Y[i] + ((1+(i*dt)^3)/3+(2/(1+(i*dt)))*Y[i])*dt + (1+(i*dt)^3)/3*sqrt(dt)*W[i]
  }
  return(Y[N+1])
}
q1 = function(X0,Y0){
  Xt = function(t){
    theta = c(.2,-.5,2/3)
    unlist(rerun(10^4,sdeGeneratorX(T=t,x0 = X0,theta = theta)))
  }
  Yt = function(t){
    unlist(rerun(10^3,sdeGeneratorY(T=t,y0=Y0)))
}
  e1 = mean(nthroot(Xt(2),3))
  e2 = mean(Yt(3))
  Xt2 = Xt(2)
  indicator = ifelse(Xt2>1,1,0)
  e3 = mean(Xt2*Yt(2)*indicator)
  p1 = mean(Yt(2)>5)
  return(rbind(p1,e1,e2,e3))
}
table1 = q1(1,3/4)
rownames(table1) = c("p1","e1","e2","e3")
colnames(table1) = "Values"
knitr::kable(table1)
Values
p1 0.9800000
e1 0.6298357
e2 26.5365293
e3 3.9850264

Question 2

sdeGenerator2 = function(N=1000,T=1,x0,theta=c(),func.type){
  t = seq(0,T,1/N)
  dt = 1/N
  X = c()
  #set.seed(123)
  W = rnorm(N)
  set.seed(1)
  Z = rnorm(N)
  if(func.type =="sde")
  for (i in 1:N) {
    X[1] =x0
    X[i+1] = X[i] + theta[1]*X[i]*dt + theta[2]*sqrt(dt)*Z[i]+theta[3]*sqrt(dt)*W[i]
    Xn = X[N+1]
  }else if(func.type =="normal"){
      Xn = mean(exp(theta[1]*T+theta[2]*W+theta[3]*Z))
    }
  return(Xn)
}
q2 = function(X0){
  Xt = function(t){
    theta = c(.25,1/3,-3/4)
     unlist(rerun(10^4,sdeGenerator2(T=t,x0 = X0,theta = theta,func.type = "sde")))
  }
  Yt = function(t){
    theta = c(-0.08,1/3,3/4)
    sdeGenerator2(T=t,x0=X0,theta=theta,func.type = "normal")
  }
  e1 = mean(nthroot((1+Xt(3)),3))
  e2 = mean(Xt(1)*Yt(1))
  return(rbind(e1,e2))
}
table2 = q2(1)
rownames(table2) = c("e1","e2")
colnames(table2) = "Values"
knitr::kable(table2)
Values
e1 1.375254
e2 2.116718

Question 3

Here the parameters using are \(S_0 = 100, T = 1, X = 100, r = 0.04, \sigma = 0.2\):

mc.call = function(S0,T,X,r,sigma,dt=0.004){
  dt = 0.004
  n = T/dt
  W = sqrt(T)*rnorm(n)
  plus=c()
  minus=c()
  payoff.plus=c()
  payoff.minus=c()
  for (i in 1:n) {
    plus[i]= S0*exp(sigma*W[i]+(r-sigma^2/2)*T)-X
    if(plus[i]>0){
      payoff.plus[i] = exp(-r*T)*(plus[i]) 
    }else{
      payoff.plus[i] = 0
    }
  }
  
  for (i in 1:n) {
    minus[i]= S0*exp(sigma*(-W[i])+(r-sigma^2/2)*T)-X
    if(minus[i]>0){
      payoff.minus[i] = exp(-r*T)*(minus[i]) 
    }else{
      payoff.minus[i] = 0
    }
  }
  
  payoff = (payoff.plus+payoff.minus)/2
  rvc = mean(payoff)
  return(rvc)
}
cat("The call option price computed by Monte Carlos is ",mc.call(100,1,100,0.04,0.2))
## The call option price computed by Monte Carlos is  10.22302
  1. Here the parameters using are \(S_0 = 100, T = 1, X = 100, r = 0.04, \sigma = 0.2\)
compute_cdf = function(x){
  if(x >= 0){
    N = 1-1/2*(1+0.0498673470*x+0.0211410061*x^2+0.0032776263*x^3+0.0000380036*x^4+0.0000488906*x^5+0.0000053830*x^6)^(-16)
  }else{
    N=1-compute_cdf(-x)#(1-1/2*(1-0.0498673470*x+0.0211410061*x^2-0.0032776263*x^3+0.0000380036*x^4-0.0000488906*x^5+0.0000053830*x^6)^(-16))
  }
  return(N)
}
bs.call = function(S0,T,K,r,sigma,dt=0.004){
  d1 = (log(S0/K) + (r + 0.5*sigma^2)*T)/(sigma*sqrt(T))
  d2 = d1 -sigma*sqrt(T)
  call = S0*compute_cdf(d1)-K*exp(-r*T)*compute_cdf(d2)
  return(call)
}
cat("The call option price computed by Black Scholes is ",bs.call(100,1,100,0.04,0.2))
## The call option price computed by Black Scholes is  9.925075
  1. Here central difference is used in estimating greeks for greater accruacy
getDelta = function(S0,T,X,r,sigma,dt=0.04,epsilon=0.0001){
  (bs.call(S0+epsilon,T,X,r,sigma,dt=0.004)-bs.call(S0-epsilon,T,X,r,sigma,dt=0.004))/(2*epsilon)
}
getGamma = function(S0,T,X,r,sigma,dt=0.04,epsilon=0.0001){
  (getDelta(S0+epsilon,T,X,r,sigma,dt=0.004)-getDelta(S0-epsilon,T,X,r,sigma,dt=0.004))/(2*epsilon)
}
getVega = function(S0,T,X,r,sigma,dt=0.04,epsilon=0.0001){
  (bs.call(S0,T,X,r,sigma+epsilon,dt=0.004)-bs.call(S0,T,X,r,sigma-epsilon,dt=0.004))/(2*epsilon)
}
getTheta  = function(S0,T,X,r,sigma,dt=0.04,epsilon=0.0001){
  -(bs.call(S0,T+epsilon,X,r,sigma,dt=0.004)-bs.call(S0,T-epsilon,X,r,sigma,dt=0.004))/(2*epsilon)
}
getGreeks = function(S0,T,X,r,sigma,dt=0.04){
  S0 =seq(15,25,1)
  Delta = c()
  Gamma = c()
  Vega = c()
  Theta = c()
  for (i in 1:length(S0)) {
    Delta[i] = getDelta(S0[i],T,X,r,sigma,dt=0.04)
    Gamma[i] = getGamma(S0[i],T,X,r,sigma,dt=0.04)
    Vega[i] = getVega(S0[i],T,X,r,sigma,dt=0.04)
    Theta[i] = getTheta(S0[i],T,X,r,sigma,dt=0.04)
  }
  return(rbind(Delta,Gamma,Vega,Theta))
  
}
S0 =seq(15,25,1)
result = getGreeks(S0,0.5,20,0.04,0.25)
par(mfrow=c(2,2))
plot(x=S0,y=result[1,],type="l",col=4,lwd=2,main = "Delta of Call",xlab="S0",ylab="Delta")
plot(x=S0,y=result[2,],type="l",col=4,lwd=2,main = "Gamma of Call",xlab="S0",ylab="Gamma")
plot(x=S0,y=result[3,],type="l",col=4,lwd=2,main = "Vega of Call",xlab="S0",ylab="Vega")
plot(x=S0,y=result[4,],type="l",col=4,lwd=2,main = "Theta of Call",xlab="S0",ylab="Theta")

Question 4

Here is this calcualtion, \(K=48\) was introduced to get a numerical solution

TwoFactorGenerator = function(rho,r,S0,V0,sigma,alpha,beta){
  N=1000
  dt = 1/N
  mu = c(0,0)
  vcov = matrix(c(1,rho,rho,1),ncol=2)
  bino = mvrnorm(N,mu,vcov)
  V.ref = numeric()
  S.ref = c()
  V.par = c()
  S.par = c()
  V.full = c()
  S.full = c()
  for (i in 1:N) {
    S.ref[1] = S.par[1] = S.full[1] = S0
    V.ref[1] = V.par[1] = V.full[1] = V0
    V.ref[i+1] = abs(V.ref[i]) + alpha*(beta-abs(V.ref[i]))*dt + sigma*sqrt(abs(V.ref[i]))*sqrt(dt)*bino[i,2]
    S.ref[i+1] = S.ref[i] + r*S.ref[i]*dt + sqrt(abs(V.ref[i]))*sqrt(dt)*bino[i,1]
    V.par[i+1] = V.par[i] + alpha*(beta-V.par[i])*dt + sigma*sqrt(pmax(0,max(V.par[i])))*sqrt(dt)*bino[i,2]
    S.par[i+1] = S.par[i] + r*S.par[i]*dt + sqrt(pmax(0,max(V.par[i])))*sqrt(dt)*bino[i,1]
    V.full[i+1] = V.full[i] + alpha*(beta-pmax(0,max(V.full[i])))*dt +sigma*sqrt(pmax(0,max(V.full[i])))*sqrt(dt)*bino[i,2]
    S.full[i+1] = S.full[i] + r*S.full[i]*dt + sqrt(pmax(0,max(V.full[i])))*sqrt(dt)*bino[i,1]
  }
  S0_1path = cbind(S.ref,S.par,S.full)[N+1,]
  return(S0_1path)
}
TwoFactorModel = function(rho,r,S0,V0,sigma,alpha,beta){
  S0_multis = rerun(10^3,TwoFactorGenerator(rho,r,S0,V0,sigma,alpha,beta))
  unlist.S0 = matrix(do.call("rbind",S0_multis),ncol=3)
  S0 = colMeans(unlist.S0)
  X=48
  call_ref = exp(-r*T)*(pmax(S0[3]-X,0)) 
  call_par = exp(-r*T)*(pmax(S0[2]-X,0))
  call_full = exp(-r*T)*(pmax(S0[1]-X,0))
  call.values = rbind(call_full,call_par,call_ref)
  return(cbind(call.values,S0))#
}
table.twofac = TwoFactorModel(-0.6,0.04,48,0.05,0.42,5.8,0.0625)
rownames(table.twofac) = c("Full Truncation","Partial Truncation","Reflection methods")
colnames(table.twofac) = c("Value of Call","Price Estimates")
knitr::kable(table.twofac)
Value of Call Price Estimates
Full Truncation 1.882206 49.95902
Partial Truncation 1.882206 49.95902
Reflection methods 1.882206 49.95902

Question 5

  1. Generate 100 2-dimensional Uniform [0,1]x[0,1] vectors by using any one of the algorithms for random number generation
LGM = function(N){
#initialize random vector
random=numeric()
#initialize uniform vector
randu.manual=numeric() 
# LGM parameters
m=2^31-1 
a=7^5
b=0
# start generation
for (i in 1:N) {
  random[1]=1
  random[i+1]=(a*random[i]+b)%% (m) # x_n+1 = (a*x_n+b) mod m
  
}
randu.manual = random/m
return(randu.manual)
}
sequence = rerun(2,LGM(100))
TwoDimensionalUniform = matrix(do.call("rbind",sequence),ncol=2)
  1. Generate 100 points of the 2-dimensional Halton sequences, using bases 2 and 7.
HaltonGenerator = function(n,m){
    output <- c()
    for (i in 1:n) {
        f = 1
        r = 0
        while (i > 0) {
            f = f/m
            r = r + f * (i%%m)
            i = floor(i/m)
        }
        output <- c(output,r)
    }
    return(output)
}
halton1 =HaltonGenerator(100,2)
halton2 =HaltonGenerator(100,7)
TwoDimensionalHalton27 = cbind(halton1,halton2)

c)Generate 100 points of the 2-dimensional Halton sequences, using bases 2 and 4. (Note: 4 is a nonprime number!).

halton1 =HaltonGenerator(100,2)
halton2 =HaltonGenerator(100,4)
TwoDimensionalHalton24=cbind(halton1,halton2)
par(mfrow=c(2,2))
plot(TwoDimensionalUniform,main="Uniform Distribution Pairs",col=2,xlab="Uniform 1",ylab="Uniform 2",pch=20)
plot(TwoDimensionalHalton27,main="Halton sequences, using bases 2 and 7",col=3,xlab="Base 2",ylab="Base 7",pch=20)
plot(TwoDimensionalHalton24,main="Halton sequences, using bases 2 and 4",col=4,xlab="Base 2",ylab="Base 4",pch=20)

Yes. There are differences. For base 2, the point distributions are more disperse and random, base 7 is more neat and uniform. For non-prime number 4, we can oberserve a linear relationship between dot distribution, and dots are more tend to bind together. The possible reason for this is that since 4 is not a prime number and 2 is a factor of it.

eval.I = function(b1,b2){
  n=10000
  x=HaltonGenerator(n,b1)
  y=HaltonGenerator(n,b2)
  Int = rep(0,n)
  for (i in 1:n) {
    Int[i] = exp(-x[i]*y[i])*(sin(6*pi*x[i])+nthroot(cos(2*pi*y[i]),3))
  }
  I = mean(Int)
  return(I)
}
eval.table = rbind(eval.I(2,4),eval.I(2,7),eval.I(5,7))
rownames(eval.table) =c("pair of 2 and 4","pair of 2 and 7","pair of 5 and 7")
colnames(eval.table) = "I"
knitr::kable(eval.table)
I
pair of 2 and 4 -0.0048839
pair of 2 and 7 0.0261144
pair of 5 and 7 0.0261637