##################
## TD ENNONCE
##################


timeVector = c(0,1,2,3)
simuVector = c(1:8)
simulationActifs <- matrix(c(1,1,1,1,1,1,1,1,1.09, 1.16, 1.22, 0.93, 1.11, 0.76, 0.92, 0.88, 1.08, 1.26, 1.07, 0.97, 1.56, 0.77, 0.84, 1.22, 1.34, 1.54, 1.03, 0.92, 1.52, 0.9, 1.01, 1.34), nrow = 8, dimnames = list(simuVector, timeVector))
simulationActifs <- as.data.frame(simulationActifs)
R = 6/100
stopping_rule = rep(3,8)
Strike = 1.1


#########################
# START THE FOR LOOP HERE
#########################
asset_at_stopping_rule = unlist(lapply(simuVector, function(x) simulationActifs[x,4]))
# Use the correct discount factor
Y = pmax(Strike-asset_at_stopping_rule,0) * exp(-R)
X = simulationActifs[[3]]

payoff_current_timeIndex = pmax(Strike-simulationActifs[[3]],0)

FILTER_ITM_PATH = (payoff_current_timeIndex > 0)

ITM_X = X[FILTER_ITM_PATH]
ITM_Y = Y[FILTER_ITM_PATH]

ORDER = order(ITM_X)

plot(ITM_X[ORDER],ITM_Y[ORDER], type="l")
model <- lm(Y[FILTER_ITM_PATH] ~ X[FILTER_ITM_PATH] + I(X[FILTER_ITM_PATH]^2))
print(model)
Yinterp = model$coefficients[1] + model$coefficients[2] * ITM_X + model$coefficients[3] * ITM_X^2
lines(ITM_X[ORDER], Yinterp[ORDER], col="red", lty = 2)

A = (payoff_current_timeIndex[FILTER_ITM_PATH] > Yinterp)
stopping_rule[FILTER_ITM_PATH][A] = 2
#########################
# START THE FOR LOOP HERE
#########################
# COMPUTE THE AMERICAN PUT
#########################
	
	
##################
## TD CORRIGE
##################


timeVector = c(0,1,2,3)
simuVector = c(1:8)
simulationActifs <- matrix(c(1,1,1,1,1,1,1,1,1.09, 1.16, 1.22, 0.93, 1.11, 0.76, 0.92, 0.88, 1.08, 1.26, 1.07, 0.97, 1.56, 0.77, 0.84, 1.22, 1.34, 1.54, 1.03, 0.92, 1.52, 0.9, 1.01, 1.34), nrow = 8, dimnames = list(simuVector, timeVector))
simulationActifs <- as.data.frame(simulationActifs)
R = 6/100
Strike = 1.1


stopping_rule = rep(3,8)
for(i in 2:1) {

	asset_at_stopping_rule = unlist(lapply(simuVector, function(x) simulationActifs[x,stopping_rule[x]+1]))
	# Use the correct discount factor
	Y = pmax(Strike-asset_at_stopping_rule,0) * exp(-R*(stopping_rule - i))
	X = simulationActifs[[i+1]]
	
	payoff_current_timeIndex = pmax(Strike-simulationActifs[[i+1]],0)

	FILTER_ITM_PATH = (payoff_current_timeIndex > 0)

	ITM_X = X[FILTER_ITM_PATH]
	ITM_Y = Y[FILTER_ITM_PATH]

	ORDER = order(ITM_X)

	plot(ITM_X[ORDER],ITM_Y[ORDER], type="l")
	model <- lm(Y[FILTER_ITM_PATH] ~ X[FILTER_ITM_PATH] + I(X[FILTER_ITM_PATH]^2))
	print(model)
	Yinterp = model$coefficients[1] + model$coefficients[2] * ITM_X + model$coefficients[3] * ITM_X^2
	lines(ITM_X[ORDER], Yinterp[ORDER], col="red", lty = 2)

	A = (payoff_current_timeIndex[FILTER_ITM_PATH] > Yinterp)
	stopping_rule[FILTER_ITM_PATH][A] = i
}

asset_at_stopping_rule = unlist(lapply(simuVector, function(x) simulationActifs[x,stopping_rule[x]+1]))
payoff_at_zero = pmax(Strike-asset_at_stopping_rule,0) * exp(-R*(stopping_rule))
american_call_price = mean(payoff_at_zero)
american_call_price




##################
## CAS REEL DE PRICING
##################

R = 6/100
Strike = 1.1

#Brownian generation
Nobs = 50
#nb d'axe à diffuser
NbDiffu = 10000
#exprimé en année dans le TD. Correspond à la maturité de la modélisation
Maturity = 3
m1 <- matrix(runif(Nobs*NbDiffu), Nobs, NbDiffu)
W = as.data.frame(apply(sqrt(Maturity/Nobs) * qnorm(m1), 2, cumsum))
if(Nobs == 1) { W = t(W) }
Wplot = rbind(rep(0, NbDiffu ), W)
xScale = Maturity*(0:Nobs)/Nobs
matplot(xScale, Wplot  , type = rep("l",NbDiffu), lty = 1, lwd = 1, col = 1:NbDiffu ,ylim  = c(-3,3))
#Spot diffusion
vol = 11/100
#Simulation de cours de bourse
S0 = 1
S = S0 * exp (  (R - vol^2/2) * xScale  + vol * Wplot )
matplot(xScale, S , type = rep("l",NbDiffu), lty = 1, lwd =1 , col =0:NbDiffu )
simulationActifs <- as.data.frame(t(S))

Nobs = Nobs+1

deltaT = Maturity/(Nobs-1)

payoff_call = function(x, Strike){
	return(pmax(x - Strike,0))
}

payoff_put = function(x, Strike){
	return(pmax(Strike - x,0))
}

chosen_payoff = payoff_put
stopping_rule = rep(Nobs-1,NbDiffu)
for(i in (Nobs-2):1) {

	asset_at_stopping_rule = unlist(lapply((1:NbDiffu), function(x) simulationActifs[x,stopping_rule[x]+1]))
	
	Y = chosen_payoff(asset_at_stopping_rule, Strike) * exp(-R*(stopping_rule - i)*deltaT)
	X = simulationActifs[[i+1]]
	
	payoff_current_timeIndex = chosen_payoff(X,Strike)

	FILTER_ITM_PATH = (payoff_current_timeIndex > 0)

	ITM_X = X[FILTER_ITM_PATH]
	ITM_Y = Y[FILTER_ITM_PATH]

	ORDER = order(ITM_X)

	# plot(ITM_X[ORDER],ITM_Y[ORDER], pch=4)
	model <- lm(Y[FILTER_ITM_PATH] ~ X[FILTER_ITM_PATH] + I(X[FILTER_ITM_PATH]^2))
	# print(model)
	Yinterp = model$coefficients[1] + model$coefficients[2] * ITM_X + model$coefficients[3] * ITM_X^2
	# lines(ITM_X[ORDER], Yinterp[ORDER], col="red", lwd = 2)
	# lines(ITM_X[ORDER], chosen_payoff(ITM_X[ORDER], Strike), col="purple", lwd = 2)

	Ai = (payoff_current_timeIndex[FILTER_ITM_PATH] > Yinterp)
	stopping_rule[FILTER_ITM_PATH][Ai] = i
	# readline(prompt="Press [enter] to continue")

}

plot(sort(stopping_rule))

asset_at_stopping_rule = unlist(lapply((1:NbDiffu), function(x) simulationActifs[x,stopping_rule[x]+1]))
payoff_at_zero = chosen_payoff(asset_at_stopping_rule,Strike) * exp(-R*(stopping_rule)*deltaT)
american_option_price = mean(payoff_at_zero)
american_option_price

d1 = (log(S0 / Strike) + (R + vol^2/2) * Maturity) / (vol * sqrt(Maturity))
d2 = (log(S0 / Strike) + (R - vol^2/2) * Maturity) / (vol * sqrt(Maturity))
CallBlackScholes = S0 * pnorm( d1 ) - Strike * exp(-R*Maturity) * pnorm( d2 )
PutBlackScholes = Strike * exp(-R*Maturity) * pnorm( -d2 ) - S0 * pnorm( -d1 ) 
CallBlackScholes
PutBlackScholes
CallBlackScholes - PutBlackScholes