# Uebungen vom 3. Blatt

# P1 
# First, we need the exponential random walk.  We will define 
# it for given parameters u,d,r,p,T, start.  

CRR_walk <- function(u,d,p,T, start){
	A = array(0, T+1)   #R counts from 1 to T+1, not from 0 to T.   
	A[1] = start
	for(i in 1:T){
		if(runif(1) < p)
		     A[i+1] = u*A[i] 
		else 
		     A[i+1] = d*A[i]
	}
	A
}

# Next, we implement, for a given path, K, r, the payoff evaluation 
# of a European option.  

European_CRR <- function(path, r, K){
	     T <- length(path)
	     end_value <- max(0, path[T] - K)*exp(-r*T)
}

# And a Monte Carlo direct simulation.  We return all results so that 
# we can also look at empricial variance.  

Eur_MC <- function(u,d,r,p,T,start,K,N){
	result = 1:N
	for(i in 1:N){
	      	path <- CRR_walk(u,d,p,T,start)
		result[i] <- European_CRR(path,r,K)
	}
	result
}


# Finally, the up-and-out:  

Upnout_CRR <- function(u,d,r,p,T,start,K,L){
	path  <- CRR_walk(u,d,p,T,start)
	if(max(path) <= L) 
	     payoff <- European_CRR(path,r,K)
	else 
	     payoff <- 0
	payoff
}


#And the Monte Carlo direct simulation 

Upnout_MC <- function(u,d,r,p,T,start,K,L,N){
	result = 1:N
	for(i in 1:N){
		result[i] =  Upnout_CRR(u,d,r,p,T,start,K,L)
	}
	result
}

# It's tedious to write this formula, so we shove it into a subroutine:  

martingal <- function(u,d,r){
	p <- (1+r-d)/(u-d)
	p
}

# Now, let's test this with some values.  

start=1
u = 1.01 
d = 0.99
r = 6/100000
T = 1000
K = 0.8
p = martingal(u,d,r)


#x = Eur_MC(u,d,r,p,T,start, K,N)
#y = Upnout_MC(u,d,r,p,T,start,K,L,N)

# Now we need to locate L such that the mean payoff is roughly 
# half the payoff of the european call.  Two problems arise:  
# How to search for L, and how to scale N, the number of experiments.  
# L will be searched after in binary search; for this we want to 
# be relatively sure that we are on the correct side with our MC guess.  
# By rough heuristic, we assume that our error is not more than 3* the 
# standard deviation, and so we will just take, as confidence interval, 
# mean plus/minus std err.  

MC_stderr <- function(data){
	result <- 3 * sqrt(var(data))/sqrt(length(data))
	result
}

# Next we will, with some care, estimate the average 
# return of the MC option.  We need a new input, the relative 
# accuracy err we wish to establish.  First we estimate mean and  
# variance by a preliminary run with N=400.  Then, we set 
# N as 6 times 4 var/(eps * mean)^2 (our goal is that  
# sqrt(var)/sqrt(N)mean <= eps/2  ), and test whether this suffices.  
# If not, increase by 2, and redo.  


Avreturn_MC <- function(u,d,r,p,T, start,K,err){
	N=400
	curr_run = Eur_MC(u,d,r,p,T,start, K,N)
	estvar = var(curr_run)
	estmean = mean(curr_run)
	if(estmean == 0){ #might happen if K too large 
	mean(curr_run)
	}
	else{
	N = ceiling(10 * 4 * estvar /(err*err*estmean*estmean))+2
	curr_run = Eur_MC(u,d,r,p,T,start, K,N)
	}
	while((mean(curr_run) != 0) && 
             (MC_stderr(curr_run)/mean(curr_run) > err/2)){
		N = ceiling(1.5*N)
	curr_run = Eur_MC(u,d,r,p,T,start, K,N)
	}
	mean(curr_run)
}

# Just for testing, the same stuff for the up and out.  Note that 
# run time scales essentially the same.  This is mainly due to the 
# fact that we always simulate the whole path - which would be 
# quite suboptimal for the european option.  

Avretuo_MC <- function(u,d,r,p,T, start,K,L, err){
	N=400
	curr_run = Upnout_MC(u,d,r,p,T,start, K,L,N)
	estvar = var(curr_run)
	estmean = mean(curr_run)
	if(estmean == 0){ #might happen if K too large or L too small 
	mean(curr_run)
	}
	else{
	N = ceiling(10 * 4 * estvar /(err*err*estmean*estmean))+2

	curr_run = Upnout_MC(u,d,r,p,T,start, K,L,N)
	}
	while((mean(curr_run) != 0) && 
              (MC_stderr(curr_run)/mean(curr_run) > err/2)){
		N = ceiling(1.5*N)
	curr_run = Upnout_MC(u,d,r,p,T,start, K,L,N)
	}
	mean(curr_run)
}


# Now we can proceed as follows:  We first estimate the average 
# payoff for european.  Then we will start with N=100 a binary 
# search.  Each time our confidence interval contains half the 
# european price, we shall increase N 6-fold.  We do this until 
# the relative accuracy of L is down to the same niveau as the 
# relative accuracy for the european payoff.  


calibrate_upnout <- function(u,d,r,p,T,start,K,err){
	europ = Avreturn_MC(u,d,r,p,T,start,K,err) #av. payoff, treated 
	goal = europ/2					 #as exact 
	if(goal == 0){L = 0}
	else{
  	 N = 100
	 L = 10*K #some start value;  
	 stepsize = L  
	 nextsim = Upnout_MC(u,d,r,p,T,start,K,L,N)   #first trial

	 while((mean(nextsim) - MC_stderr(nextsim) < goal) && 
		(mean(nextsim) + MC_stderr(nextsim) > goal)){
		N = 6*N
	 nextsim = Upnout_MC(u,d,r,p,T,start,K,L,N)
	 }			#refine (increase N) if necessary
	 
	 while(stepsize/L > err/2){
		stepsize = stepsize/2
	  	if( (mean(nextsim) - MC_stderr(nextsim)) > goal){
		    L = L - stepsize  #The binary search itself
		}
		else {L = L + stepsize}
  		nextsim = Upnout_MC(u,d,r,p,T,start,K,L,N)

		while(mean(nextsim) - MC_stderr(nextsim) < goal && 
		        mean(nextsim) + MC_stderr(nextsim) > goal){
			N = 6*N
			nextsim = Upnout_MC(u,d,r,p,T,start,K,L,N)
		}  #Again, some refinement if necessary  
	 }
	} 
	alarm()
	L
}

# Now, on to some fun.  Let's compile for a bunch (say, twenty) of 
# different strike prices upper bounds, plot them, and so on.  
# Warning:  This can take some hours to run if you set average 
# accuracy to something useful like 0.05.  

uno_vs_euro <- function(u,d,r,p,T,start,range,err){
	    results = array(0, length(range)) 
	    for(j in 1:length(range)){
	    	  k = range[j]
	    	  results[j] = calibrate_upnout(u,d,r,p,T,start,k,err)
	    }
	    results
}


# Our next project looks for the number of fixpoint-free permutations 
# in S_{1000}.  

n = 1000 

isfixfree <- function(perm) {
	result = 1
	for(i in 1:length(perm)){ 
		if( perm[i] == i) {result = 0} 
	result
	}
}

 # Now we recycle our method for generating random permutations (see 
 # ueb1.R):   

draw <- function(k,n){
	result <- array(0, k)
	X <- 1:n      #This is the drawing range 
	for(i in 1:k){
		m <- discuni(1,n-i+1)
	result[i] <- X[m]   #Now X[m] was drawn and should be kicked out.  
	X <- X[X != X[m]]   #awkward kludge since R has no proper slicing
	}	
	result 
}

# P(fixp.free) = #fixp.free / 1000!.  So we approximate the LHS 
# with Monte Carlo:  

M = 100000 

fixp_MC <- function(n,M){
	sum = 0; 
	for(i in 1:M){
	perm = draw(n,n)
	sum = sum + isfixfree(perm) #count one up if fixp.fre
	}
	sum = sum/M
	sum
}

#Now we would need to multiply by an ungodly number, 1000!.  We skip 
#this.  Instead, we plot, for some values of n, the extimated relative 
#number of fixpoint-free permutations:  

fixfreeplot <- function(n,M){
	x = array(0, n)
	for(i in 1:n){
	x[i] <- fixp_MC(50*i,M)
	}
	x
}

# Last we shall simulate equidistribution on a simplex, namely 
# on x_i >= 0, x_1 + .. + x_d <= 1.  

# The straightforward idea is direct rejection:  

reject <- function(d){
	x <- runif(d)
	while(sum(x) > 1){x <- runif(d)} #repeat until in simplex
	x
}

# Well, this is.. slow.  A little thinking reveals:  If x is uniform 
# on the simplex, then, given x_1, .., x_i, x_{i+1} is uniform on 
# 0..(1 - (x_1 + .. x_i)).  So the following is a much faster method:  

scale <- function(d){
	x <- runif(d)
	for(i in 2:d){
	x[i] = x[i] *(1 -  sum(x[1:(i-1)]))
	}
	x
}

# At first sight, it seems strange that almost all results have 
# sum > 0.99; this is, however, to be expected, since the overwhelming
# part of the volume lies in the 'fat' part of the simplex; in the 
# extremal case of an infinite-dimensional simplex, all mass lies on 
# the border sum = 1.