# Skripty v tomto souboru pochazi z knihy (s drobnymi upravami): 
# Stefano M. Iacus: Option Pricing and Estimation of Financial Models with R. John Wiley & Sons, 2011. 



AmericanPutExplicit <- function (Smin = 0, Smax, T = 1, N = 10, M = 10, K, r = 0.05, sigma = 0.01) {
	Dt <- T / N
	DS <- (Smax - Smin) / M
	t <- seq (0, T, by = Dt)
	S <- seq (Smin, Smax, by = DS)
	A <- function (j) {
		(-0.5 * r * j * Dt + 0.5 * sigma^2 * j^2 * Dt) / (1 + r * Dt)
	}
	B <- function (j) {
		(1 - sigma^2 * j^2 * Dt)/(1 + r * Dt)
	}
	C <- function (j) {
		(0.5 * r * j * Dt + 0.5 * sigma^2 * j^2 *	Dt) / (1 + r * Dt)
	}
	P <- matrix (, M + 1, N + 1)
	colnames (P) <- round (t, 2)
	rownames (P) <- round (rev(S), 2)
	P[M + 1,] <- K
	P[1, ] <- 0
	P[, N + 1] <- sapply (rev (S), function (x) {max(K - x, 0)})
	optTime <- matrix (FALSE, M + 1, N + 1)
	optTime[M+1, ] <- TRUE
	optTime[which (P[, N+1] > 0), N+1] <- TRUE
	for (i in (N - 1):0) {
		for (j in 1:(M - 1)) {
			J <- M + 1 - j
			I <- i + 1
			P[J, I] <- A (j) * P[J+1, I+1] + B (j) * P[J, I+1] + C(j) * P[J-1, I+1]
			if (P[J, I] < P[J, N+1]) {optTime[J, I] <- TRUE}
		}
	}
	colnames (optTime) <- colnames (P)
	rownames (optTime) <- rownames (P)
	ans <- list (P = P, t = t, S = S, optTime = optTime, N = N, M = M)
	class (ans) <- "AmericanPut"
	return (invisible (ans))
}

plot.AmericanPut <- function(obj) {
	plot (range (obj$S), range (obj$t), type = "n", axes = F, xlab = "S", ylab = "t")
	axis (1, obj$S, obj$S)
	axis (2, obj$t, obj$t)
	abline (h = obj$t, v = obj$S, col = "darkgray", lty = "dotted")
	for (i in 0:obj$N) {
		for (j in 0:obj$M) {
			J <- obj$M + 1 - j
			I <- i + 1
			cl <- "red"
			if (obj$optTime[J, I]) {cl <- "blue"}
			text (obj$S[j+1], obj$t[i+1], round (obj$P[J, I], 1), cex = 0.75, col = cl)
		}
	}
	DS <- mean (obj$S[1:2])
	y <- as.numeric (apply (obj$optTime, 2, function (x) {which(x)[1]}))
	lines (obj$S[obj$M+2-y] + DS, obj$t, lty = 2)
}



# ==============================================================================



put <- AmericanPutExp(Smin=0, Smax=60, T=1, N=10, M=10, K=30, r=0.05, sigma=0.01)

put$P

plot (put)



# ==============================================================================



AmericanPutImplicit <- function (Smin = 0, Smax, T = 1, N = 10, M = 10, K, r = 0.05, sigma = 0.01) {
	Dt <- T / N
	DS <- (Smax - Smin) / M
	t <- seq (0, T, by = Dt)
	S <- seq (Smin, Smax, by = DS)
	A <- function (j) {
		0.5 * r * j * Dt - 0.5 * sigma^2 * j^2 * Dt
	}
	B <- function (j) {
		1 + sigma^2 * j^2 * Dt + r * Dt
	}
	C <- function (j) {
		-0.5 * r * j * Dt - 0.5 * sigma^2 * j^2 *	Dt
	}
	a <- sapply (0:M, A)
	b <- sapply (0:M, B)
	c <- sapply (0:M, C)
	P <- matrix (, M + 1, N + 1)
	colnames (P) <- round (t, 2)
	rownames (P) <- round (rev(S), 2)
	P[M+1, ] <- K
	P[1, ] <- 0
	P[, N+1] <- sapply (rev(S), function (x) {max (K - x, 0)})
	AA <- matrix(0, M - 1, M - 1)
	for (j in 1:(M - 1)) {
		if (j > 1) {
			AA[j, j-1] <- A (j)
		}
		if (j < M) {
			AA[j, j] <- B (j)
		}
		if (j < M - 1) {
			AA[j, j+1] <- C (j)
		}
	}
	optTime <- matrix (FALSE, M+1, N+1)
	for (i in (N - 1):0) {
		I <- i + 1
		bb <- P[M:2, I+1]
		bb[1] <- bb[1] - A(1) * P[M+1-0, I+1]
		bb[M-1] <- bb[M-1] - C (M-1) * P[M+1-M, I+1]
		P[M:2, I] <- solve (AA, bb)
		idx <- which (P[, I] < P[, N+1])
		P[idx, I] <- P[idx, N+1]
		optTime[idx, I] <- TRUE
	}
	optTime[M+1, ] <- TRUE
	optTime[which (P[, N+1] > 0), N+1] <- TRUE
	colnames (optTime) <- colnames (P)
	rownames (optTime) <- rownames (P)
	ans <- list (P = P, t = t, S = S, optTime = optTime, N = N, M= M)
	class (ans) <- "AmericanPut"
	return (invisible (ans))
}



# ==============================================================================



put2 <- AmericanPutImp(Smax = 60, sigma = 0.4, K = 30)

put2$P

plot (put2)



# ==============================================================================



library (fOptions)

opt <- BAWAmericanApproxOption("p", S=36, X=30, Time=1, r=0.05, b=0.05, sigma=0.4)
opt@price

?BSAmericanApproxOption



# ==============================================================================