Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Calibration of financial models
########## 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)
########## 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)
########## 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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.