0
$\begingroup$

I am working with a Heston model discretization through truncation, given by the following code:

(for (i in 1:Nsteps){
X<-log(S)
X<-X+(R-0.5*pmax(V,0))*dt+sqrt(pmax(V,0))*dW1
S<-exp(X)
V<-V+kappa*(theta-pmax(V,0))*dt+sigma*sqrt(pmax(V,0))*dW2

}))

I am trying to simulate several paths let's say 10 or 100 of the stock under Hestons model above. I have used the function replicate in the following fashion:

P_A_cap_B<-replicate(10,(for (i in 1:Nsteps){
X<-log(S)
X<-X+(R-0.5*pmax(V,0))*dt+sqrt(pmax(V,0))*dW1
S<-exp(X)
V<-V+kappa*(theta-pmax(V,0))*dt+sigma*sqrt(pmax(V,0))*dW2

}))

However I am getting null values in each collum. To summarize it is not working with the for loop

If I do:

replicate(10,S)

All the values across the columns are equal, which is clearly what I did not want.

Question:

How to simulate 10 different paths of the stock?

Thanks in Advance!

Appendix(Full code):

 BlackScholesFormula  <- function (spot,timetomat,strike,r, q=0, sigma, opttype=1, greektype=1)
{ 

d1<-(log(spot/strike)+ ((r-q)+0.5*sigma^2)*timetomat)/(sigma*sqrt(timetomat))
d2<-d1-sigma*sqrt(timetomat)

if (opttype==1 && greektype==1) result<-spot*exp(-q*timetomat)*pnorm(d1)-strike*exp(-r*timetomat)*pnorm(d2)

if (opttype==2 && greektype==1) result<-spot*exp(-q*timetomat)*pnorm(d1)-strike*exp(-r*timetomat)*pnorm(d2)-spot*exp(-q*timetomat)+strike*exp(-r*timetomat)

if (opttype==4 && greektype==1) result<-(spot^2)*exp((r+sigma^2)*timetomat)
if (opttype==4 && greektype==2) result<-2*spot*exp((r+sigma^2)*timetomat)

if (opttype==1 && greektype==2) result<-exp(-q*timetomat)*pnorm(d1)


if (opttype==2 && greektype==2) result<-exp(-q*timetomat)*(pnorm(d1)-1)

if (greektype==3) result<-exp(-q*timetomat)*dnorm(d1)/(spot*sigma*sqrt(timetomat))

if (greektype==4) result<-exp(-q*timetomat)*spot*dnorm(d1)*sqrt(timetomat)

BlackScholesFormula<-result

}

BlackScholesImpVol  <- function (obsprice,spot,timetomat,strike,r, q=0, opttype=1)

  { difference<- function(sigBS, obsprice,spot,timetomat,strike,r,q,opttype)
      {BlackScholesFormula (spot,timetomat,strike,r,q,sigBS, opttype,1)-obsprice
     }

    uniroot(difference, c(10^-6,10),obsprice=obsprice,spot=spot,timetomat=timetomat,strike=strike,r=r,q=q,opttype=opttype)$root

  }

BS.Fourier <- function(spot,timetoexp,strike,r,divyield,sigma,greek=1){

  X<-log(spot/strike)+(r-divyield)*timetoexp
  
  integrand<-function(k){
    integrand<-Re(exp( (-1i*k+0.5)*X- 0.5*(k^2+0.25)*sigma^2*timetoexp)/(k^2+0.25))
  } 
  dummy<-integrate(integrand,lower=-Inf,upper=Inf)$value
  
  BS.Fourier<-exp(-divyield*timetoexp)*spot-strike*exp(-r*timetoexp)*dummy/(2*pi)
}

Heston.Fourier <- function(spot,timetoexp,strike,r,divyield,V,theta,kappa,epsilon,rho,greek=1){

X<-log(spot/strike)+(r-divyield)*timetoexp
kappahat<-kappa-0.5*rho*epsilon
xiDummy<-kappahat^2+0.25*epsilon^2

integrand<-function(k){
    xi<-sqrt(k^2*epsilon^2*(1-rho^2)+2i*k*epsilon*rho*kappahat+xiDummy)
    Psi.P<--(1i*k*rho*epsilon+kappahat)+xi
    Psi.M<-(1i*k*rho*epsilon+kappahat)+xi
    alpha<--kappa*theta*(Psi.P*timetoexp + 2*log((Psi.M + Psi.P*exp(-xi*timetoexp))/(2*xi)))/epsilon^2
    beta<--(1-exp(-xi*timetoexp))/(Psi.M + Psi.P*exp(-xi*timetoexp))
    numerator<-exp((-1i*k+0.5)*X+alpha+(k^2+0.25)*beta*V)

    if(greek==1) dummy<-Re(numerator/(k^2+0.25))
    if(greek==2) dummy<-Re((0.5-1i*k)*numerator/(spot*(k^2+0.25)))
    if(greek==3) dummy<--Re(numerator/spot^2)
    if(greek==4) dummy<-Re(numerator*beta)
    
    integrand<-dummy
  } 


  dummy<-integrate(integrand,lower=-100,upper=100)$value

  if (greek==1) dummy<-exp(-divyield*timetoexp)*spot-strike*exp(-r*timetoexp)*dummy/(2*pi)

  if(greek==2) dummy<-exp(-divyield*timetoexp)-strike*exp(-r*timetoexp)*dummy/(2*pi)

  if(greek==3) dummy<--strike*exp(-r*timetoexp)*dummy/(2*pi)
  
  if(greek==4) dummy<--strike*exp(-r*timetoexp)*dummy/(2*pi)

  Heston.Fourier<-dummy
}

Andreasen.Fourier <- function(spot,timetoexp,strike,Z,lambda,beta,epsilon){

X<-log(spot/strike)

integrand<-function(k){
    neweps<-lambda*epsilon
    xi<-sqrt(k^2*neweps^2+beta^2+0.25*neweps^2)
    Psi.P<--beta+xi
    Psi.M<-beta+xi
    A<--beta*(Psi.P*timetoexp + 2*log((Psi.M + Psi.P*exp(-xi*timetoexp))/(2*xi)))/(epsilon^2)
    B<-(1-exp(-xi*timetoexp))/(Psi.M + Psi.P*exp(-xi*timetoexp))
    integrand<-Re(exp( (-1i*k+0.5)*X+A-(k^2+0.25)*B*lambda^2*Z)/(k^2+0.25))
  } 

  dummy<-integrate(integrand,lower=-Inf,upper=Inf)$value

 Andreasen.Fourier<-spot-strike*dummy/(2*pi)
  
}


timetoexp<-1.0

S0<-100
R<-0.02
V0<-0.15^2
kappa<-2
theta<-0.2^2
sigma<-1.0
rho<--0.5

Params<-paste("S0=", S0, ", sqrt(V0)=",sqrt(V0),", r=",R, ", T=", timetoexp)
ModelParams<-paste("kappa =", kappa, ", theta =", theta, ", rho =", rho, ", sigma =", sigma)


Nsim<-10^4
NstepsPerYear<-1*252

Nsteps<-round(timetoexp*NstepsPerYear)
dt<-timetoexp/Nsteps

S<-rep(S0,Nsim); V<-rep(V0,Nsim)

AbsAtZero<-FALSE
if (AbsAtZero) title<-"Heston: Euler-S + |V| @0"
MaxAtZero<-FALSE
if (MaxAtZero) title<-"Heston: Euler-S + V^+ @0"
FullTruncation<-TRUE
# From page 6 in https://www.dropbox.com/s/nw7uzmf8k0imq0t/LeifHestonWP.pdf?dl=0
# "The scheme that appears to produce the smallest discretization bias
if (FullTruncation) title<-"Heston\n Euler-ln(S) + full truncation-V"
Andersen_1<-
if (Andersen_1) title<-"Heston: Euler-S + |V| @0"


ToFile<-FALSE; GraphFile="FullTruncation.png" 
IVspace<-FALSE

RunTime<-system.time(
  for (i in 1:Nsteps){
    dW1<-sqrt(dt)*rnorm(Nsim,0,1)
    dW2<-rho*dW1+sqrt(1-rho^2)*sqrt(dt)*rnorm(Nsim,0,1)
    if (AbsAtZero){ 
      S<-S+R*S*dt+sqrt(V)*S*dW1
      V<-abs(V + kappa*(theta-V)*dt + sigma*sqrt(V)*dW2)
    }
    if (MaxAtZero){ 
      S<-S+R*S*dt+sqrt(V)*S*dW1
      V<-pmax(V + kappa*(theta-V)*dt + sigma*sqrt(V)*dW2,0)
    }
    
    if (FullTruncation){
      X<-log(S)
      X<-X+(R-0.5*pmax(V,0))*dt+sqrt(pmax(V,0))*dW1
      S<-exp(X)
      V<-V+kappa*(theta-pmax(V,0))*dt+sigma*sqrt(pmax(V,0))*dW2
    }
    
    if(Andersen_1){
      S<-S+R*S*dt+sqrt(V)*S*dW1    
      alpha_h<-0.5
      gamma_1<-sqrt((1/dt)*log(1+((0.5*sigma*sigma)*(1/kappa)*V^(2*alpha_h)*(1-exp(-2*kappa*dt)))/((exp(-kappa*dt)*V+(1-exp(-kappa*dt))*theta)^2)))
      V<-(exp(-kappa*dt)*V+(1-exp(-kappa*dt))*theta)*exp(-0.5*(gamma_1*gamma_1)+gamma_1*dW2)
    }
    
  }  
)[3]

ExperimentParams<-paste("#steps/year =",NstepsPerYear, ", #paths =", Nsim, ", runtime (seconds) = ", round(RunTime,2))


StdErrSimCall<-TrueCall<-SimCall<-Strikes<-S0+(-10:10)*2
ConfBandIVSimCall<-IVSimCall<-IVTrueCall<-SimCall

if(ToFile) png(GraphFile,width=14,height=14,units='cm',res=300)

for (i in 1:length(SimCall)) {
  SimCall[i]<-exp(-R*timetoexp)*mean(pmax(S-Strikes[i],0))
  StdErrSimCall[i]<-sd(exp(-R*timetoexp)*(pmax(S-Strikes[i],0)))/sqrt(Nsim)
  if (IVspace) IVSimCall[i]<-BlackScholesImpVol(SimCall[i],S0,timetoexp,Strikes[i],R, q=0, opttype=1)
  if (IVspace) ConfBandIVSimCall[i]<-BlackScholesImpVol(SimCall[i]+1.96*StdErrSimCall[i],S0,timetoexp,Strikes[i],R, q=0, opttype=1)
  TrueCall[i]<-Heston.Fourier(S0,timetoexp,Strikes[i],R,0,V0,theta,kappa,sigma,rho,greek=1)
  if (IVspace) IVTrueCall[i]<-BlackScholesImpVol(TrueCall[i],S0,timetoexp,Strikes[i],R, q=0, opttype=1)
}

if(!IVspace){
  dummy<-c(min(SimCall,TrueCall),max(SimCall,TrueCall))
  plot(Strikes,SimCall,type='b',ylim=dummy,col='blue',ylab="Call price",main=title,xlab="Strike")
  points(Strikes,TrueCall,type='l')
  points(Strikes,SimCall+1.96*StdErrSimCall,type='l',lty=2,col='blue')
  points(Strikes,SimCall-1.96*StdErrSimCall,type='l',lty=2,col='blue')
  
  text(min(Strikes),min(dummy)+0.10*(max(dummy)-min(dummy)),"Black: Closef-form Heston ala Lipton",adj=0,cex=0.5)
  text(min(Strikes),min(dummy)+0.05*(max(dummy)-min(dummy)),"Blue: Sim' est' (o) + 95% conf' band' (---)", col="blue",adj=0,cex=0.5)
  
  text(max(Strikes),max(dummy),Params,adj=1, cex=0.5)
  text(max(Strikes),0.95*max(dummy),ModelParams,adj=1, cex=0.5)
  text(max(Strikes),0.90*max(dummy),ExperimentParams,adj=1, cex=0.5)
  text(max(Strikes),0.85*max(dummy),".../Dropbox/FinKont2/LiveAtLectures_HestonSimulation.R",adj=1, cex=0.5)
}

if(IVspace){
  dummy<-c(min(2*IVSimCall-ConfBandIVSimCall,IVTrueCall),max(IVTrueCall,ConfBandIVSimCall) )
  plot(Strikes,IVSimCall,type='b',col='blue',ylim=dummy, ylab="Implied volatility",main=title,xlab="Strike")
  points(Strikes,IVTrueCall,type='l')
  points(Strikes,ConfBandIVSimCall,type='l',lty=2,col='blue')
  points(Strikes,2*IVSimCall-ConfBandIVSimCall,type='l',lty=2,col='blue')
  
  text(min(Strikes),min(dummy)+0.10*(max(dummy)-min(dummy)),"Black: Closef-form Heston ala Lipton",adj=0,cex=0.5)
  text(min(Strikes),min(dummy)+0.05*(max(dummy)-min(dummy)),"Blue: Sim' est' (o) + 95% conf' band' (---)", col="blue",adj=0,cex=0.5)
  
  text(max(Strikes),max(dummy),Params,adj=1, cex=0.5)
  text(max(Strikes),0.95*max(dummy),ModelParams,adj=1, cex=0.5)
  text(max(Strikes),0.90*max(dummy),ExperimentParams,adj=1, cex=0.5)
  text(max(Strikes),0.85*max(dummy),".../Dropbox/FinKont2/LiveAtLectures_HestonSimulation.R",adj=1, cex=0.5)
}

if (ToFile)
  
# Control Variates
  
  
   S_vc<-rep(S0,Nsim)
  
for (i in 1:Nsteps){
    dW1<-sqrt(dt)*rnorm(Nsim,0,1)
    dW2<-rho*dW1+sqrt(1-rho^2)*sqrt(dt)*rnorm(Nsim,0,1)
    S_vc<-S_vc+R*S_vc*dt+sqrt(V0)*S_vc*dW1}


print(S_vc)

#Consider strike


strike_v0<-75

#equation 6)P_A_CV=P_A_cap+(P_G-P_G_cap)

#Parameters
#Estimating P_A_cap
P_A_cap_vec<-rep(0,length(S))
for (i in 1:length(S)){P_A_cap_vec[i]<-(max((S[i]-strike_v0),0))}
print(P_A_cap_vec) 
P_A_cap<-mean(P_A_cap_vec)
print(P_A_cap)

#Estimating P_G
P_G<-BlackScholesFormula(S0,timetoexp,strike_v0,R,0,V0,1,1)
print(P_G)
#Estimating P_G_cap

P_G_cap_vec<-rep(0,length(S_vc))
for(i in 1:length(S_vc)){
P_G_cap_vec<-(max(S_vc[i]-strike_v0,0))}
P_G_cap<-mean(P_G_cap_vec)
print(P_G_cap)

#Equation 6) Control Variate
P_A_CV<-P_A_cap+(P_G-P_G_cap)

#Estiamtor of the Heston model Option Pricing
print(P_A_CV)

#equation 7)P_A_CV=P_A_cap+Beta*(P_G-P_G_cap), using the parameter beta to minimize variance

#Generate random S and S_vc paths and compute P_A_cap and P_G_cap for each path.

set.seed(190)

P_A_cap_B<-replicate(10,{for (i in 1:Nsteps){
X<-log(S)
X<-X+(R-0.5*pmax(V,0))*dt+sqrt(pmax(V,0))*dW1
S<-exp(X)
V<-V+kappa*(theta-pmax(V,0))*dt+sigma*sqrt(pmax(V,0))*dW2}

})
print(P_A_cap_B)
Improve this question
$\endgroup$
4
  • $\begingroup$ Where are dW1 and dW2 coming from? $\endgroup$
    – Bob Jansen ♦
    Mar 21, 2021 at 19:18
  • $\begingroup$ @BobJansen Jansen dW1<-sqrt(dt)*rnorm(Nsim,0,1) dW2<-rho*dW1+sqrt(1-rho^2)*sqrt(dt)*rnorm(Nsim,0,1) $\endgroup$
    – Pedro Gomes
    Mar 21, 2021 at 19:32
  • $\begingroup$ Please make your code reproducible. It’s hard to help if important details are missing. $\endgroup$
    – Bob Jansen ♦
    Mar 21, 2021 at 19:43
  • $\begingroup$ @BobJansen I have just uploaded the full code in my answer. $\endgroup$
    – Pedro Gomes
    Mar 21, 2021 at 19:48

1 Answer 1

Reset to default
3
$\begingroup$

The replicate function works best when you fully define your discretization scheme within a function. Then you can simply replicate the function-call x amount of times. Also, try and keep code duplication to a minimum and improve your general syntax. This will help you and your peers that might need to review and/or change your code in the future.

Nevertheless, one way to solve your problem is to define the discretization scheme for the Heston model within a function:

HestonSimTruncation <- function(S0, V0, kappa, theta, rho, epsilon, r, Nsteps, T){

  V <- numeric() 

  S <- numeric()

  dt <- T/Nsteps 

  dw1 <- sqrt(dt) * rnorm(Nsteps) #Generation of Brownian increments for dw1

  dw2 <- rho*dw1+sqrt((1-rho^2)*dt)*rnorm(Nsteps) #Generation of Brownian increments for dw2

  S[1] <- log(S0) 

  V[1] <- V0 

  for(i in 1:Nsteps){

    V[i+1] <- V[i] + kappa*(theta-V[i])*dt + epsilon*sqrt(V[i]) * dw1[i]

    if( V[i+1] < 0 ){ V[i+1] <- 0 } #Absorption barrier, similar to full truncation. 

    S[i+1] = S[i] + (r-0.5 * V[i]) * dt + sqrt(V[i]) * dw2[i]

  } 

  return(exp(S))

}

Then you can input your parameters and simply do:

test2 <- replicate(10, HestonSimTruncation(S0, V0, kappa, theta, rho, epsilon, r, Nsteps, T))

Which will give you a matrix with a dimension of $\mathbb{R}^{Nsteps \times Npaths}$, where $Npaths = 10$. Below, I have provided a graphical illustration of the first 5 simulations using the above function, hestonsims

where I have used the parameters $S_0 = 100$, $V_0 = 0.2^2$, $\kappa = 3.8$, $\theta = 0.3095$, $\rho = -0.78$, $\epsilon = 0.9288$, $Nsteps = 1000$, $T=1$ (I found these parameters on the internet). I hope this help.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.