# R uses per default the Mersenne twister, which is rather good.  
# A very simple test, just for demonstration, is Kolmogorov-Smirnov:  

x = runif(100000)
ks.test(x, "punif")

# A casino game.  With probability p, you win your investment each round, 
# with probability 1-p, you lose your investment.  You start with S, you 
# end iff S <= 0 or S >= K (if S << 0, you might have some real problems).  
# We first study the strategy to invest q*(K-S) + r for q, r >= 0.  

# First the casino round; we return the difference in money.  

round <- function(x, p){
	if(runif(1) < p) result = x 
	else result = -x
	result
}

# Now the first strategy; return is 1 for success and 0 for failure 

casino1 <- function(p,S,K,q,r){
	while(S > 0 && S < K){
	S = S + round( min(q*(K-S) + r, S), p)
	}
	if(S<=0) {result = 0} else {result = 1}
	result
}

# We will make a direct MC simulation, for S=100, K=200 and p=0.45.  

S <- 100 
K <- 200
p <- 0.4

strat1 <- function(q,r,N){

	succ = 0
	for(n in 1:N){
		succ = succ + casino1(p,S,K,q,r)
	}
	succ/N
}

# A simple insurance model.  The payoff process is a Poisson process.  

# Following the exercise, we model this with waiting times.  
# In each step, we look for the next jump; if it is > T, we 
# survived, otherwise we look at the capital at the jump time.  
# Since, obviously, the company's maximal loss occurs at the 
# jump time, this tells us if we survive until the next jump time.  

# First a single survival decision:  

survival <- function(lambda, T, rho, K){
	t <- 0
	survived <- 1
	cval <- K + 1
	wait <- 0 
	while( ((t < T) && (survived == 1))){
		cval = cval - 1 + rho * wait 
		if(cval < 0) {survived <- 0; break}
		wait <- rexp(1, rate=lambda)
		t <- t + wait
	}
	survived
}

# Now a direct MC simulation 

MCsurvival <- function(lambda, T, rho, K, N){
	nosurv <- 0
	for(i in 1:N){
	nosurv <- nosurv + survival(lambda, T, rho, K)
	}
	nosurv/N
}


# How to create random permutations?  The simplest idea is to 
# create a random element from {1,..,n}^n, then check whether it 
# is a permutation.  This, however, is insanely costly.  Our first 
# method demonstrates this:  

discuni <- function(k, n){
	result <- ceiling(runif(k,min=0, max=n))
	result
}

naiveperm <- function(n){
	isperm <- 0	
	while(isperm == 0){
		x = discuni(n,n)
		isperm <- 1
		y <- sort(x)
		for(i in 2:n){
			if(y[i] == y[i-1])
			isperm <- 0
		}		
	}
	x
}


# A more sophisticated idea is to check on each new number if 
# it is already used:  

perm2 <- function(n){
	x = array(0, n)
	for(i in 1:n){
		isnew <- 0 
		while(isnew == 0){
			new = discuni(1,n)
			if(length(x[x==new])==0) #means:  new is not in x yet
			isnew <- 1
			}
		x[i] <- new
	}
	x
}

# But this is still slow; on average, one needs ~n^2 trials.  A better 
# idea is to implement 'draw k of n' directly.  This is done by shrinking 
# the 'drawing range' after each draw by kicking out the drawn element.  

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 
}

#Still, this is a bit clumsy, since slicing means shifting of a lot of 
#stuff in memory.  Here is a more elegant way, using ranks.  If X1,..,Xn 
#are iid continuously distributed, and if r(n) is the position of the nth 
#largest element, then r(1),..,r(n) is a random permutation, and this yields 
#the uniform distribution.  The rank is easily found as the position in the 
#order statistics, and already implemented in R.