<jupyter_start><jupyter_text>AFV code demonstrationWe demonstrate the code used in the paper *Efficient simulation of affine forward variance models*.<jupyter_code>library(repr)
source("BlackFormula.R")
source("AFVsimulation.R")
source("GammaKernel.R")
source("roughHestonAdams.R")
source("roughHestonPade.R")
source("Lewis.R")
bl <- "royalblue"
rd <- "red2"
pk <- "hotpink1"
gr <- "green4"<jupyter_output><empty_output><jupyter_text>Parameters We choose parameters similar to those found from a fit to SPX options as of May 19, 2017, the same data that was used in Roughening Heston.<jupyter_code>params0 <- list(al=0.55,lam=0,eta=0.8,rho=-0.65, H=0.05,lam=0)
xi0 <- function(s){0.025+0*s} # The forward variance curve<jupyter_output><empty_output><jupyter_text>Simulation using the RSQE and HQE schemes The option "all" returns a list with many variables of interest.<jupyter_code>system.time(rsqe.sim <- RSQE.sim(params0, xi0)(paths=100000, steps=100, expiries=1,output="all"))
system.time(hqe.sim.100 <- HQE.sim(params0, xi0)(paths=100000, steps=100, expiries=1,output="all"))
system.time(hqe.sim.200 <- HQE.sim(params0, xi0)(paths=100000, steps=200, expiries=1,output="all"))
names(hqe.sim.100)<jupyter_output><empty_output><jupyter_text>- v is terminal variance $v_T$- x is terminal log-spot $X_T$- y is terminal value of $Y_T = \int_0^T\,\sqrt{v_s}\,dW_s$.- w is terminal value of quadratic variation $w_T = \int_0^T\,{v_s}\,ds$.  ParallelizationThe code can be run in parallel with selected variables as output.For example, here with $v_T$ as output.<jupyter_code>library(foreach)
library(doParallel)
paths <- 1000000 # Note 1 million paths here!
steps <- 100
t0<-proc.time()

# Number of iterations
iters<- max(1,floor(paths/1000))

# Setup parallel backend to use all processors
(cl.num <- detectCores()) # Display number of processors on your machine
cl<-makeCluster(cl.num)
registerDoParallel(cl)

# Loop
ls <- foreach(icount(iters)) %dopar% {
        HQE.sim(params0, xi0)(paths=1000, steps=steps, expiries=1,output="v")
        }

stopCluster(cl)
    hqe.sim.v <- do.call(cbind, ls) #Bind all of the submatrices into one big matrix

print(proc.time()- t0)<jupyter_output><empty_output><jupyter_text>Figure 7.6: Histograms of $v_T$<jupyter_code>options(repr.plot.width=14,repr.plot.height=7,repr.plot.res=150)
vv <- hqe.sim.v
vvg <- vv[vv > 1e-9] # Restrict sample to values greater than tiny
par(mfrow=c(1,2))
hist(vv,breaks=200,xlab=expression(v[T]),main="",col=bl, border=bl, cex.lab=1.5)
par(new=F)
hist(vvg,breaks=200,xlab=expression(v[T]),main="",col=bl,border=bl, cex.lab=1.5)
par(mfrow=c(1,1))<jupyter_output><empty_output><jupyter_text>Histograms of terminal stock prices<jupyter_code>spots.100 <- exp(hqe.sim.100$x)
spots.200 <- exp(hqe.sim.200$x)
hist(spots.100,breaks=200,xlab=expression(S[T]),main="",col=pk, border=pk, cex.lab=1.5)<jupyter_output><empty_output><jupyter_text>Note the fat negative tail in the return distribution.  Draw the 1 year smile with the above parameters<jupyter_code>smile.100 <- function(k){ivS(spots.100, T=1, mean(spots.100)*exp(k))}
smile.200 <- function(k){ivS(spots.200, T=1, mean(spots.200)*exp(k))}<jupyter_output><empty_output><jupyter_text>We select a vector of log-strikes.<jupyter_code>kk <- seq(-.4,.4,.01)
smile.HQE.100 <- smile.100(kk)
smile.HQE.200 <- smile.200(kk)
options(repr.plot.width=10,repr.plot.height=7,repr.plot.res=150)
plot(kk,smile.HQE.100,col=rd,lwd=2,type="l",
     xlab="Log-strike k", ylab = "Implied vol.",cex.lab=1.5)
lines(kk,smile.HQE.200,col=bl,lwd=2,lty=2)
legend("topright",c("HQE 100 steps","HQE 200 steps"), cex=1.5, inset=.05, 
       lty=c(1,2),col=c(rd,bl), lwd=2)<jupyter_output><empty_output><jupyter_text>Richardson extrapolated smile<jupyter_code>smile.HQE.Richardson <- 2*smile.200(kk)-smile.100(kk)
plot(kk,smile.HQE.100,col=rd,lwd=2,type="l",
     xlab="Log-strike k", ylab = "Implied vol.",cex.lab=1.5)
lines(kk,smile.HQE.200,col=bl,lwd=2,lty=2)
lines(kk,smile.HQE.Richardson,col=gr,lwd=2,lty=4)
legend("topright",c("HQE 100 steps","HQE 200 steps",
    "HQE Richardson"), cex=1.5, inset=.05, 
       lty=c(1,2,4),col=c(rd,bl,gr), lwd=2)<jupyter_output><empty_output><jupyter_text>Parameter conversionThe QE simulation and the Adams/Padé approximations use different formulations of the rough Heston model.In the QE case,$$d\xi_t(u) = \eta\,\sqrt{2 H}\,(u-t)^{\alpha-1}.$$In the Adams/Padé case,$$d\xi_t(u) = \frac{\nu}{\Gamma(\alpha)}\,(u-t)^{\alpha-1}.$$Thus$$\nu = \sqrt{2 H}\, \Gamma(\alpha)\,\eta.$$<jupyter_code>(params0$nu <- params0$eta*sqrt(2*params0$al-1)*gamma(params0$al))<jupyter_output><empty_output><jupyter_text>Comparison with Adams smile - The Adams scheme code is orignally due to <NAME>, <NAME>, <NAME>, and <NAME>, University of Bologna (2020).- The code for computing the cutoff in the Lewis formula is based on code originally by <NAME>, École Polytechnique Paris (2017). Note that the Adams smile takes time to compute (two minutes on my machine).<jupyter_code>system.time(vol.Adams.kk.1y.200 <-  impliedVolRoughHeston(params0, xi0, nSteps=200)(kk,1))
plot(kk,smile.HQE.Richardson,col=rd,lwd=2,type="l",
     xlab="Log-strike k", ylab = "Implied vol.",cex.lab=1.5)
lines(kk,vol.Adams.kk.1y.200,col=bl,lwd=2,lty=2)
legend("topright",c("HQE Richardon","Adams"), cex=1.5, inset=.05, 
       lty=c(1,2),col=c(rd,bl), lwd=2)<jupyter_output><empty_output><jupyter_text>Padé approximation<jupyter_code>phi.Pade <- phiRoughHestonDhApprox(params0, xi0, dh.approx= d.h.Pade33, n=500)
system.time(vol.Pade.kk.1y <-  impvol.phi(phi.Pade)(kk,1))
plot(kk,smile.HQE.Richardson,col=rd,lwd=2,type="l",
     xlab="Log-strike k", ylab = "Implied vol.",cex.lab=1.5)
lines(kk,vol.Adams.kk.1y.200,col=bl,lwd=2,lty=2)
lines(kk,vol.Pade.kk.1y,col=gr,lwd=2,lty=4)
legend("topright",c("HQ<NAME>","Adams","Padé"), cex=1.5, inset=.05, 
       lty=c(1,2,2),col=c(rd,bl,gr), lwd=2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1807412421729565, 'fraction_numerical': 0.053985445929937384, 'mean_word_length': 6.415307402760352, 'pattern_counts': {'":': 0, '<': 86, '<?xml version=': 0, '>': 56, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '581058', 'license': 'MIT', 'n_tokens_mistral': 2672, 'n_tokens_neox': 2400, 'n_words': 426, 'path': 'AFV simulation demonstration.ipynb', 'repo_name': 'jgatheral/AFVsimulationCode', 'size': '5956'}