{{ message }}

Instantly share code, notes, and snippets.

WilsonMongwe/Fitting a Jump Diffusion.r

Last active Sep 30, 2017
Calibration of financial models
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 ########## 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 file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 ########## 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 file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 ########## 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")