Calibration of financial models

Fitting a Jump Diffusion.r

########## This code calculates non parametric estimates of the jump-difusion parameters #######
simulateJump=function(mu_,ss_,lambda_,mu2_,sigma_,TotalTime,delta)
{
  Sn=0
  times <- c(0)
  while(Sn <= TotalTime)
  {
    n <- length(times)
    u <- runif(1)
    expon <- -log(u)/lambda_
    Sn <- times[n]+expon
    times <- c(times, Sn)
  }
  
  times=times[-length(times)]   
  t=seq(from=0,to=TotalTime,by=delta)
  #the last time is beyond TotalTime, so i delete from the vector times
  indicator=seq(from=0,to=0,length=length(t))  # if there are jumps or not between two times
  jumpSize=seq(from=0,to=0,length=length(t))   # stores the size of the jump
  
 
  for(m in 2:length(times))
  {
    for(k in 2: length(t))
    {
      if( t[k-1]<=times[m] && times[m]<=t[k])
      {
        indicator[k]=1
        jumpSize[k]=rnorm(1,mean=mu2_,sd=sigma_)
      }
    }
  }
  
  stock=seq(from=0,to=0,length=length(t))
  stock[1]=10
  
  for(i in 2:length(t))
  {
    stock[i]=stock[i-1]+(mu_)*delta+ss_*sqrt(delta)*rnorm(1)+jumpSize[i]*indicator[i]
  }
  
  final=cbind(t,stock,jumpSize,indicator)
  return(final)
}

#### define the moment functions

M1 <- function(a, xdifference,timediff,xcurrentValue, h)
{
  kernel.weights <- K( (xcurrentValue-a)/h )
  xdiff <-xdifference/(timediff/delta)
  sum(kernel.weights*xdiff)/(delta*sum(kernel.weights)) 
}


M2 <- function(a, xdifference,timediff,xcurrentValue, h)
{
  kernel.weights <- K( (xcurrentValue-a)/h )
  xdiff <- ( xdifference/(timediff/delta) )^2
  sum(kernel.weights*xdiff)/( delta*sum(kernel.weights) ) 
}


M4<- function(a, xdifference,timediff,xcurrentValue, h)
{
  kernel.weights <- K( (xcurrentValue-a)/h )
  xdiff <-( xdifference/(timediff/delta) )^4
  sum(kernel.weights*xdiff)/(delta*sum(kernel.weights)) 
}

M6 <- function(a, xdifference,timediff,xcurrentValue, h)
{
  kernel.weights <- K( (xcurrentValue-a)/h )
  xdiff <-( xdifference/(timediff/delta) )^6
  sum(kernel.weights*xdiff)/(delta*sum(kernel.weights)) 
}


delta=1/(8*252) # stock exchages are open on average for about 8 hours every businessday, and there are 252 business days.
mu_=0.0 # drift
ss_=0.15 # jump diffusion diffusive vol
lambda_=25# jump intensity
mu2_=0 # mean jump size 
sigma_=0.03 # jump size volatility
TotalTime=20# years. In this case 20 years

dataMax=simulateJump(mu_=mu_,ss_=ss_,lambda_=lambda_,mu2_=mu2_,sigma_=sigma_,TotalTime=TotalTime,delta=delta)
data=dataMax[,1:2]

Xdiff = diff(data[,2])
Xcurrent = data[1:length(Xdiff),2]
TimeDiff = diff(data[,1])

V=cbind(Xdiff,TimeDiff,Xcurrent)

#### define the kernel function, in this case being the Gausian kernel

K <- function(u)
{
  exp(-0.5*u^2)/sqrt(2*pi)
}

############ Calculating the estimates for the parameters ##########
n=200
X=Xcurrent
H=sort(X)   #values at which we are evaluating the institaneous

# The bandwiths for the moments
h1= 3
h2= 2
h4= 1
h6= 0.9

#first.est <-sapply(H, M1, xdifference=V[,1],timediff=V[,2], xcurrentValue=V[,3], h=h1)
second.est <-sapply(H, M2, xdifference=V[,1],timediff=V[,2], xcurrentValue=V[,3], h=h2)
fourth.est <-sapply(H, M4, xdifference=V[,1],timediff=V[,2], xcurrentValue=V[,3], h=h4)
sixth.est <-sapply(H, M6, xdifference=V[,1],timediff=V[,2], xcurrentValue=V[,3], h=h6)

############ Calculating the bootstrap estimates  ############

first.boot.est=c()
second.boot.est=c()
forth.boot.est=c()
sixth.boot.est=c()

for(i in 1:n)
{
  sampleRows<- sample(1:length(V[,1]), size=length(V[,1]), replace=TRUE)
  NewV=V[sampleRows,]
  
  #first.boot.est<- cbind(first.boot.est, sapply(H, M1,xdifference=NewV[,1],timediff=NewV[,2],xcurrentValue=NewV[,3],h=h1)) 
  second.boot.est<- cbind(second.boot.est, sapply(H, M2,xdifference=NewV[,1],timediff=NewV[,2],xcurrentValue=NewV[,3],h=h2)) 
  forth.boot.est<- cbind(forth.boot.est, sapply(H, M4,xdifference=NewV[,1],timediff=NewV[,2],xcurrentValue=NewV[,3],h=h4)) 
  sixth.boot.est<- cbind(sixth.boot.est, sapply(H, M6,xdifference=NewV[,1],timediff=NewV[,2],xcurrentValue=NewV[,3],h=h6))
  NewV=c()
  print(paste("Iteration Number :",i))
}



####### plots of the estimated parameters   ########

### plot for drift

plot(H, first.est, type="l", ylim=c(-0.5,0.5), cex.lab=1.5, lwd=2, ylab=expression(mu), xlab="Log share price")
for(i in 1:n) {  lines(H,first.boot.est[,i],col="grey") }

### confidence interval for drift

lower.first=rep(0,length(H))
upper.first=rep(0,length(H))
for(k in 1:length(H))
{ lower.first[k]=quantile(first.boot.est[k,],prob=c(0.1))
upper.first[k]=quantile(first.boot.est[k,],prob=c(0.9))
}
lines(H,upper.first,lty="dashed",lwd=2,col="blue")
lines(H,first.est,lty="dashed",lwd=2,col="green")
lines(H,lower.first,lty="dashed",lwd=2,col="red")
legend("topright", legend=c("90% boostrap bound", "Estimate","10% boostrap bound"),
       col=c("red", "green","blue"), lty=1:2, cex=0.8, text.font=4)

###### Estimatong the jump variance, intensity function and sigma ########

## estimates of  jump variance, intensity function and sigma using the original data ##

var_y=mean(sixth.est/(5*fourth.est))

sqrt(var_y) # This is the estimate of the jump size vol

lam=fourth.est/(3*var_y^2) # Estimate of jump intensity

sig=(second.est-lam*var_y )

## bootstrap estimates of jump variance, intensity function and sigma ##

lamvec=c()
variancevec=c()
sigmavec=c()

for( i in 1:n)
{
  var_y.boot=mean(sixth.boot.est[,i]/(5*forth.boot.est[,i]))
  lam.boot=forth.boot.est[,i]/(3*(var_y.boot)^2)
  sig.boot= (second.boot.est[,i]-lam.boot*var_y.boot )
  
  variancevec=rbind(variancevec,var_y.boot)
  lamvec=cbind(lamvec,lam.boot)
  sigmavec=cbind(sigmavec,sig.boot)
  
}

##### confidence interval for the jump size standard deviation####
quantile(sqrt(variancevec),probs=c(0.1,0.9))

########### plots of lamda and sigma ###########

### plot for lambda###

plot(H, lam, type="l",ylim=c(15,35),cex.lab=1.5,lwd=2, ylab=expression(lambda), xlab="Log share price")
for(i in 1:n) {  lines(H,lamvec[,i],col="grey") }

### confidence interval for lambda ###

lower.lambda=rep(0,length(H))
upper.lambda=rep(0,length(H))
for(k in 1:length(H))
{ lower.lambda[k]=quantile(lamvec[k,],prob=c(0.1))
upper.lambda[k]=quantile(lamvec[k,],prob=c(0.9))
}

lines(H,upper.lambda,lty="dashed",lwd=2,col="blue")
lines(H,lam,lty="dashed",lwd=2,col="green")
lines(H,lower.lambda,lty="dashed",lwd=2,col="red")
legend("topright", legend=c("90% boostrap bound", "Estimate","10% boostrap bound"),
       col=c("blue", "green","red"), lty=1:2, cex=0.8, text.font=4)


### plot for sigma ###
plot(H, sqrt(sig), type="l",ylim=c(0.1,0.2),lwd=2, cex.lab=1.5, ylab=expression(sigma[jump-diffusion]), xlab="Log share price")
for(i in 1:n) {  lines(H,sqrt(sigmavec[,i]),col="grey") }

### confidence interval for sigma ###

lower.sig=rep(0,length(H))
upper.sig=rep(0,length(H))
for(k in 1:length(H))
{ lower.sig[k]=quantile(sigmavec[k,],prob=c(0.1))
upper.sig[k]=quantile(sigmavec[k,],prob=c(0.9))
}

lines(H,sqrt(upper.sig),lty="dashed",lwd=2,col="blue")
lines(H,sqrt(sig),lty="dashed",lwd=2,col="green")
lines(H,sqrt(lower.sig),lty="dashed",lwd=2,col="red")
legend("topright", legend=c("90% boostrap bound", "Estimate","10% boostrap bound"),
       col=c("blue", "green","red"), lty=1:2, cex=0.8, text.font=4)

Fitting a Standard Diffsuion.r

########## This code calculates non parametric estimates of the standard diffusion parameters #######

simulateGBM=function(mu_,ss_,TotalTime,delta)
{
  time=seq(from=0,to=TotalTime,by=delta)
  stock=seq(from=0,to=0,length=length(time))
  
  stock[1]=10
  
  for(i in 2:length(time))
  {
    stock[i]=stock[i-1]+(mu_)*delta+ss_*sqrt(delta)*rnorm(1)
  }
  
  final=cbind(time,stock)
  return(final)
}

#### define the sample moment functions
M1 <- function(a, xdifference,timediff,xcurrentValue, h)
{
  kernel.weights <- K( (xcurrentValue-a)/h )
  xdiff <-xdifference/(timediff/delta)
  sum(kernel.weights*xdiff)/(delta*sum(kernel.weights)) 
}


M2 <- function(a, xdifference,timediff,xcurrentValue, h)
{
  kernel.weights <- K( (xcurrentValue-a)/h )
  xdiff <- ( xdifference/(timediff/delta) )^2
  sum(kernel.weights*xdiff)/( delta*sum(kernel.weights) ) 
}

### generate the data #####
delta=1/(8*252) # This is hourly data per day for 252 business days
mu_= 0.055 # drift
ss_= 0.305 # diffusive vol
TotalTime= 20 # years. In this case 20 years

data=simulateGBM(mu_=mu_,ss_=ss_,TotalTime=TotalTime,delta=delta)

Xdiff = diff(data[,2])
Xcurrent = data[1:length(Xdiff),2]
TimeDiff = diff(data[,1])

V=cbind(Xdiff,TimeDiff,Xcurrent)

#### define the kernel fucntion, in this case being the Gausian kernel

K <- function(u)
{
  exp(-0.5*u^2)/sqrt(2*pi)
}


############ Calculating the estimates for the parameters ##########
n=200  # number of boostrap samples
X=Xcurrent
H=sort(X)   #values at which we are evaluating the institaneous

# The bandwiths for the moments
h1 = 3
h2 = 2

first.est <-sapply(H, M1, xdifference=V[,1],timediff=V[,2], xcurrentValue=V[,3], h=h1)
second.est <-sapply(H, M2, xdifference=V[,1],timediff=V[,2], xcurrentValue=V[,3], h=h2)

############ Calculating the bootstrap estimates  ############

first.boot.est=c()
second.boot.est=c()
forth.boot.est=c()
sixth.boot.est=c()

for(i in 1:n)
{
  sampleRows<- sample(1:length(V[,1]), size=length(V[,1]), replace=TRUE)
  NewV=V[sampleRows,]

  first.boot.est<- cbind(first.boot.est, sapply(H, M1,xdifference=NewV[,1],timediff=NewV[,2],xcurrentValue=NewV[,3],h=h1)) 
  second.boot.est<- cbind(second.boot.est, sapply(H, M2,xdifference=NewV[,1],timediff=NewV[,2],xcurrentValue=NewV[,3],h=h2)) 
  
  NewV=c()
  
  print(paste("Iteration Number :",i))
}

####### plots of the estimated parameters   ########

### plot for drift

plot(H, first.est, type="l", ylim=c(-1,1), cex.lab=1.5, lwd=2, ylab=expression(mu), xlab="Log share price")
for(i in 1:n) {  lines(H,first.boot.est[,i],col="grey") }

### confidence interval for drift

lower.first=rep(0,length(H))
upper.first=rep(0,length(H))
for(k in 1:length(H))
{ lower.first[k]=quantile(first.boot.est[k,],prob=c(0.1))
  upper.first[k]=quantile(first.boot.est[k,],prob=c(0.9))
}

lines(H,upper.first,lty="dashed",lwd=2,col="blue")
lines(H,first.est,lty="dashed",lwd=2,col="green")
lines(H,lower.first,lty="dashed",lwd=2,col="red")
legend("topright", legend=c("90% boostrap bound", "Estimate","10% boostrap bound"),
       col=c("red", "green","blue"), lty=1:2, cex=0.8, text.font=4)

### plot for the diffsuive volatility

plot(H, sqrt(second.est), cex.lab=1.5, type="l",ylim=c(0.28,0.32),lwd=2, ylab=expression(hat(sigma)[standard-diffusion]), xlab="Log share price")
for(i in 1:n) {  lines(H,sqrt(second.boot.est[,i]),col="grey") }

### confidence interval for var

lower.second=rep(0,length(H))
upper.second=rep(0,length(H))
for(k in 1:length(H))
{ lower.second[k]=quantile(sqrt(second.boot.est[k,]),prob=c(0.1))
  upper.second[k]=quantile(sqrt(second.boot.est[k,]),prob=c(0.9))
}


lines(H,upper.second,lty="dashed",lwd=2,col="blue")
lines(H,sqrt(second.est),lty="dashed",lwd=2,col="green")
lines(H,lower.second,lty="dashed",lwd=2,col="red")

legend("topright", legend=c("90% boostrap bound", "Estimate","10% boostrap bound"),
       col=c("red", "green","blue"), lty=1:2, cex=0.8, text.font=4)

Fitting the JSE Top 40.r

########## This code calculates non parametric estimates of the parameters #######

### read in data from excell file Data for model fitting
data=read.table("clipboard",header=TRUE)
attach(data)
Xdiff
Xcurrent
TimeDiff

V=cbind(Xdiff,TimeDiff,Xcurrent)


#### define the kernel fucntion, in this case being the Gausian kernel

K <- function(u)
{
  exp(-0.5*u^2)/sqrt(2*pi)
}

#### define the moment functions

M1 <- function(a, xdifference,timediff,xcurrentValue, h)
{
  kernel.weights <- K( (xcurrentValue-a)/h )
  xdiff <-xdifference/(timediff/delta)
   sum(kernel.weights*xdiff)/(delta*sum(kernel.weights)) 
}


M2 <- function(a, xdifference,timediff,xcurrentValue, h)
{
  kernel.weights <- K( (xcurrentValue-a)/h )
  xdiff <- ( xdifference/(timediff/delta) )^2
   sum(kernel.weights*xdiff)/( delta*sum(kernel.weights) ) 
}


M4<- function(a, xdifference,timediff,xcurrentValue, h)
{
  kernel.weights <- K( (xcurrentValue-a)/h )
  xdiff <-( xdifference/(timediff/delta) )^4
  sum(kernel.weights*xdiff)/(delta*sum(kernel.weights)) 
}

M6 <- function(a, xdifference,timediff,xcurrentValue, h)
{
  kernel.weights <- K( (xcurrentValue-a)/h )
  xdiff <-( xdifference/(timediff/delta) )^6
   sum(kernel.weights*xdiff)/(delta*sum(kernel.weights)) 
}


############ Calculating the estimates for the parameters ##########
n=1000

delta=15/(24*60) # delta is 15 minute
X=Xcurrent
H=sort(X)   #values at which we are evaluating the institaneous

# The bandwiths for the moments
h1=0.031
h2=0.022
h4=0.027
h6=0.027

first.est <-sapply(H, M1, xdifference=V[,1],timediff=V[,2], xcurrentValue=V[,3], h=h1)
second.est <-sapply(H, M2, xdifference=V[,1],timediff=V[,2], xcurrentValue=V[,3], h=h2)
fourth.est <-sapply(H, M4, xdifference=V[,1],timediff=V[,2], xcurrentValue=V[,3], h=h4)
sixth.est <-sapply(H, M6, xdifference=V[,1],timediff=V[,2], xcurrentValue=V[,3], h=h6)


############ Calculating the bootstrap estimates  ############

first.boot.est=c()
second.boot.est=c()
forth.boot.est=c()
sixth.boot.est=c()

for(i in 1:n)
{
  sampleRows<- sample(1:length(V[,1]), size=length(V[,1]), replace=TRUE)
  NewV=V[sampleRows,]

  first.boot.est<- cbind(first.boot.est, sapply(H, M1,xdifference=NewV[,1],timediff=NewV[,2],xcurrentValue=NewV[,3],h=h1)) 
  second.boot.est<- cbind(second.boot.est, sapply(H, M2,xdifference=NewV[,1],timediff=NewV[,2],xcurrentValue=NewV[,3],h=h2)) 
  forth.boot.est<- cbind(forth.boot.est, sapply(H, M4,xdifference=NewV[,1],timediff=NewV[,2],xcurrentValue=NewV[,3],h=h4)) 
  sixth.boot.est<- cbind(sixth.boot.est, sapply(H, M6,xdifference=NewV[,1],timediff=NewV[,2],xcurrentValue=NewV[,3],h=h6))
  NewV=c()
}



####### plots of the estimated parameters   ########


### plot for drift

plot(H, first.est, type="l", ylim=c(-0.05,0.05), cex.lab=1.5, lwd=2, ylab=expression(mu), xlab="Log of the JSE Top 40 index")
for(i in 1:n) {  lines(H,first.boot.est[,i],col="grey") }

### confidence interval for drift

lower.first=rep(0,length(H))
upper.first=rep(0,length(H))
for(k in 1:length(H))
{ lower.first[k]=quantile(first.boot.est[k,],prob=c(0.1))
  upper.first[k]=quantile(first.boot.est[k,],prob=c(0.9))
}
lines(H,lower.first,lty="dashed",lwd=2)
lines(H,upper.first,lty="dashed",lwd=2)


### plot for second moment

plot(H, sqrt(second.est), cex.lab=1.5, type="l",ylim=c(0.007,0.018),lwd=2, ylab=expression(hat(sigma)[standard-diffusion]), xlab="Log of the JSE Top 40 index")
for(i in 1:n) {  lines(H,sqrt(second.boot.est[,i]),col="grey") }

### confidence interval for drift

lower.second=rep(0,length(H))
upper.second=rep(0,length(H))
for(k in 1:length(H))
{ lower.second[k]=quantile(sqrt(second.boot.est[k,]),prob=c(0.1))
  upper.second[k]=quantile(sqrt(second.boot.est[k,]),prob=c(0.9))
}
lines(H,lower.second,lty="dashed",lwd=2)
lines(H,upper.second,lty="dashed",lwd=2)




###### Estimatong the jump variance, intensity function and sigma ########


## estimates of  jump variance, intensity function and sigma using the original data ##


var_y=mean(sixth.est/(5*fourth.est))
lam=fourth.est/(3*var_y^2)
sig=(second.est-lam*var_y )

## bootstrap estimates of jump variance, intensity function and sigma ##

lamvec=c()
variancevec=c()
sigmavec=c()

for( i in 1:n)
{
var_y.boot=mean(sixth.boot.est[,i]/(5*forth.boot.est[,i]))
lam.boot=forth.boot.est[,i]/(3*(var_y.boot)^2)
sig.boot= (second.boot.est[,i]-lam.boot*var_y.boot )

variancevec=rbind(variancevec,var_y.boot)
lamvec=cbind(lamvec,lam.boot)
sigmavec=cbind(sigmavec,sig.boot)

}

##### confidence interval for the jump size standard deviation####
quantile(sqrt(variancevec),probs=c(0.1,0.9))

########### plots of lamda and sigma ###########

### plot for lambda###

plot(H, lam, type="l",ylim=c(0.00,90),cex.lab=1.5,lwd=2, ylab=expression(lambda), xlab="Log of the JSE Top 40 index")
for(i in 1:20) {  lines(H,lamvec[,i],col="grey") }

### confidence interval for lambda ###


lower.lambda=rep(0,length(H))
upper.lambda=rep(0,length(H))
for(k in 1:length(H))
{ lower.lambda[k]=quantile(lamvec[k,],prob=c(0.1))
  upper.lambda[k]=quantile(lamvec[k,],prob=c(0.9))
}
lines(H,lower.lambda,lwd=2, lty="dashed")
lines(H,upper.lambda,lwd=2, lty="dashed")


### plot for sigma ###

plot(H, sig, type="l",ylim=c(0,0.00025),lwd=2, cex.lab=1.5, ylab=expression(sigma[jump-diffusion]), xlab="Log of the JSE Top 40 index")
for(i in 1:n) {  lines(H,sigmavec[,i],col="grey") }

### confidence interval for sigma ###

lower.sig=rep(0,length(H))
upper.sig=rep(0,length(H))
for(k in 1:length(H))
{ lower.sig[k]=quantile(sigmavec[k,],prob=c(0.1))
  upper.sig[k]=quantile(sigmavec[k,],prob=c(0.9))
}

plot(H, sqrt(sig), type="l",ylim=c(0,0.015),lwd=2,cex.lab=1.5, ylab=expression(sigma[jump-diffusion]), xlab="Log of the JSE Top 40 index")
lines(H,sqrt(lower.sig),lwd=2, lty="dashed")
lines(H,sqrt(upper.sig),lwd=2, lty="dashed")