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