library(MASS) library(SPDest) ## set the random seed (change it to obtain different simulated values :-) set.seed(1) ## set the number of observations npoz=50 ## set the true call and put prices and corresponding SPD true.call=cbind(c(10,20,30,40),c(17.5,9,2.5,2)) true.put=cbind(c(10,20,30,40),c(3,4.5,8,17.5)) true.spd=c(0.15,0.2,0.6,0.05) x=cbind(rbind(true.call,true.call),c(rep(0,4),rep(0,4))) ## calculate the nr. of replication of each strike and random order of transactions rep.sim=rmultinom(1,npoz,rep(true.spd/2,2)) x.sim=apply(x,2,rep,times=rep.sim) xord=order(x.sim[,3],x.sim[,1]) x.sim=x.sim[xord,] ## generate the uniform transaction times temp=spdest.xd(x.sim) ttime=sqrt(runif(npoz)) ## true covariance matrix for the random errors using function spdest.cvr() true.var=spdest.cvr(ttime,temp$L,temp$xd,nc=npoz-sum(x.sim[,3]),np=sum(x.sim[,3]),CL=TRUE,PT=FALSE) ## generate the correlated random errors tmp=eigen(true.var) var12=tmp$vectors%*%diag(sqrt(tmp$values))%*%t(tmp$vectors) errors=0.5*var12%*%as.matrix(rnorm(npoz)) pr=x.sim pr[,2]=pr[,2]+errors ## print the matrix used as input for the programme ## 1.column contains the strike prices ## 2.column contains the (random) option prices ## 3.column is the put-call indicator (here always 0 -- calls only) print(pr) ## run the programme using 3 variants of the algorithm a.full=spdest(price=pr ,ttime=ttime) a.crab=spdest(price=pr ,ttime=ttime,fast=TRUE,fast.type="crab",nrec=2,iter=2) a.brab=spdest(price=pr ,ttime=ttime,fast=TRUE,fast.type="brab",nrec=2,iter=2) ## plot the results par(mfrow=c(3,1)) plot(pr,col=1,xlab="strike",ylab="price",main=paste("SPD estimate, ",npoz," observations.",sep="")) lines(a.full$fit.call,lty=2) lines(true.call,lty=1) plot(pr,col=1,xlab="strike",ylab="price",main=paste("SPD estimate (CRAB), ",npoz," observations.",sep="")) lines(a.crab$fit.call,lty=2) lines(true.call,lty=1) plot(pr,col=1,xlab="strike",ylab="price",main=paste("SPD estimate (BRAB), ",npoz," observations.",sep="")) lines(a.brab$fit.call,lty=2) lines(true.call,lty=1)