#source("hedge.050207.R")

# INCLUDE THE B/S-formula
source(file="BlackScholesFormula.R")
#BlackScholesFormula(S(t),T-t,K,r, q=0, sigma, opttype=1, greektype=1)


# INITIALIZE
# ==========

S0<-100
r<-0.05
mu<-0.1
sigma<-0.2

capT<-1
strike<-105

Nhedge<-floor(sqrt(2)^(6:12)) 
Nrep<-10000

Vpf<-discoptionpayoff<-dischedgeerror<-matrix(ncol=length(Nhedge),nrow=Nrep)


# HEDGE
# =====

for(j in 1:length(Nhedge)){

St<-seq(S0, S0, length=Nrep)
dt<-capT/Nhedge[j]
initialoutlay<-BlackScholesFormula(S0,capT,strike, r,0,sigma,1,1)

Vpf[,j]<-initialoutlay

a<-BlackScholesFormula(St,capT,strike, r,0,sigma,1,2)
b<-Vpf[,j]-a*St

for(i in 2:Nhedge[j]){

    St<-St*exp((mu-0.5*sigma^2)*dt +sigma*sqrt(dt)*rnorm(Nrep))	
    Vpf[,j]<-a*St+b*exp(dt*r)    
    a<- BlackScholesFormula(St,(capT-(i-1)*dt),strike, r,0,sigma,1,2)
    b<-(Vpf[,j]-a*St)
  }
  
  St<-St*exp((mu-0.5*sigma^2)*dt +sigma*sqrt(dt)*rnorm(Nrep))
  Vpf[,j]<-a*St+b*exp(dt*r) 
  dischedgeerror[,j]<-exp(-r*capT)*(Vpf[,j]-pmax(St-strike,0))
  discoptionpayoff[,j]<-exp(-r*capT)*pmax(St-strike,0)
	
}

# SUMMARY STATS & GRAPHS
# ======================

stddevdiscerror<-1:length(Nhedge)

for(j in 1:length(Nhedge)){

stddevdiscerror[j]<-sqrt(var(dischedgeerror[,j]))
output<-paste("With ", Nhedge[j], " hedge points the std. dev. of the discounted hedge error is ", round(stddevdiscerror[j],5))

}

par(mfrow=c(1,1))


plot(Nhedge,stddevdiscerror/initialoutlay,log="xy", xlab="# hedge points", ylab="hedge error",ylim=c(0.01,0.3))
fit<-lm(log(stddevdiscerror/initialoutlay)~log(Nhedge))
points(Nhedge,exp(fit$coefficients[1])*Nhedge^fit$coefficients[2],type='l')
par(adj=0)
text(Nhedge[2],stddevdiscerror[1]/initialoutlay,paste("slope = ",round(fit$coefficients[2],3)))
par(adj=0.5)

plot(St,Vpf[,length(Nhedge)])