# =================================================================================
# 1 
# Wienerov proces

dt = 0.001
t = seq(0, 1, by = dt)
W0 = 0

generuj.Wp <- function(mu = 0, sigma = 1, dt = 0.001, W0 = 0, t = seq(0, 1, by = dt)) {
  N <- length(t)
  dW <- rnorm(N - 1, mu, sigma) * sqrt(dt)
  W <- cumsum(c(W0, dW))
  plot (t, W, type = "l", xlab = "t", ylab = "W")
  abline (h = W0, lty = 2)
  return(W)
}

W <- generuj.Wp()

# =================================================================================
# 2a
# N trajektorii WP

N <- 10  # pozadovany pocet trajektorii
mu <- 0
sigma <- 0.3

generuj.NWp <- function(N, mu = 0, sigma = 1, dt = 0.001, W0 = 0, t = seq(0, 1, by = dt)){
  M <- sapply (1:N , function (k) {
    generuj.Wp(mu = mu, sigma = sigma, dt = dt, W0 = W0, t = t)
  })
  matplot(t, M, type = "l", lty = 1, col = "#ccaaee", xlab = "t", ylab = "W", main = "Trajektorie W")
}

generuj.NWp(N, mu, sigma)

# =================================================================================
# 2b
# Geometricky WP

N <- 10  # pozadovany pocet trajektorii

generuj.NgBp <- function(N, mu = 0, sigma = 1, dt = 0.001, W0 = 100, t = seq(0, 1, by = dt)){
  generuj.gBp <- function (mu = mu, sigma = sigma, dt = dt, W0 = W0, t = t){
    n <- length (t)
    dW <- rnorm (n - 1, mu, sigma) * sqrt (dt)
    dX <- 1 + mu * dt + sigma * dW
    X <- cumprod (c (W0 , dX))
  }
  M <- sapply(1:N, function(k){
    generuj.gBp(mu = mu, sigma = sigma, dt = dt, W0 = W0, t = t)
  })
  matplot (t, M, type = "l", 
           lty = 1, xlab = "t", ylab = "geometricky Brownuv pohyb", col = "grey")
}

generuj.NgBp(N)

# =================================================================================
# 3

dt = 0.001
t = seq(0, 1, by = dt)
W0 = 0

# Dvojrozmerny Wienerov proces:
# nahodny
W1 <- generuj.Wp (t = t, dt = dt, W0 = W0)
W2 <- generuj.Wp (t = t, dt = dt, W0 = W0)
plot(W1, W2, type="l")
points(W1[length(W1)],W2[length(W2)], col="red", type="p", lwd=4) 

# Trojrozmerny Wienerov proces:
W1 <- generuj.Wp (t = t, dt = dt, W0 = W0)
W2 <- generuj.Wp (t = t, dt = dt, W0 = W0)
W3 <- generuj.Wp (t = t, dt = dt, W0 = W0)

# library(rgl)
plot3d(W1,W2,W3, type="l")
points3d(W1[length(W1)],W2[length(W2)],W3[length(W3)], col="red", lwd=150)

# =================================================================================
# 4
# Cena opcie

call.opce <- function(K, S, sigma, r, t){
  d1 <- (log(S/K) + (r + 0.5*sigma^2) * t ) / (sigma * sqrt(t))
  price <- S*pnorm(d1) - K*exp(-r*t)*pnorm(d1 - sigma*sqrt(t))
  return(price)
}

put.opce <- function(K, S, sigma, r, t){
  d1 <- (log(S/K) + (r + 0.5*sigma^2) * t ) / (sigma * sqrt(t))
  price <- -S*pnorm(-d1) + K*exp(-r*t)*pnorm(-(d1 - sigma*sqrt(t)))
  return(price)
}

# K <- 15; S <- 14; sigma <- 0.3; r <- 0.05; t <- 1
# call.opce(K, S, sigma, r, t)
# put.opce(K, S, sigma, r, t)

# =================================================================================
# 5
# Greeks

S <- 55
K <- 50
sigma <- 0.3
r <- 0
t <- 0.5

Delta <- function(S, K, t, r, sigma){
  d1 <- (log(S/K) + (r + 0.5*sigma^2) * t) / (sigma * sqrt (t))
  delta <- pnorm(d1)
  return(delta)
}

Gama <- function(S, K, t, r, sigma){
  d1 <- (log(S/K) + (r + 0.5*sigma^2) * t) / (sigma * sqrt (t))
  gama <- ( (1/sqrt(2*pi) * exp(-d1^2/2) ) / (S * sigma * sqrt(t)))
  return(gama)
}

Theta <- function(S, K, t, r, sigma){
  d1 <- (log(S/K) + (r + 0.5*sigma^2) * t) / (sigma * sqrt (t))
  theta <- - (S*sigma / (2*sqrt(t))) * (1/sqrt(2*pi) * exp(-d1^2/2)) 
           - K*r*exp(-r*t) * pnorm(d1 - sigma*sqrt(t))
  return(theta)
}

Delta(S, K, t, r, sigma)
Gama(S, K, t, r, sigma)
Theta(S, K, t, r, sigma)

#=================================================================================================
# 6
# Zavislost Delta na S

T1 <- 1
T2 <- 0.5
T3 <- 0.1
T4 <- 0.01

d1 <- rep(0,101)
d2 <- rep(0,101)
d3 <- rep(0,101)
d4 <- rep(0,101)

S<-seq(0,100,by=1)

for (i in 1: 101){
  d1[i] <- Delta(S[i], K, T1, r, sigma)
  d2[i] <- Delta(S[i], K, T2, r, sigma)
  d3[i] <- Delta(S[i], K, T3, r, sigma)
  d4[i] <- Delta(S[i], K, T4, r, sigma)
}
  
plot(0, 0, type="n", xlab = "S", ylab = "Delta", xlim = c(0,100), ylim = c(0,1))
lines(S, d1, col = 1)
lines(S, d2, col = 2)
lines(S, d3, col = 3)
lines(S, d4, col = 4)
legend("bottomright", c("T1","T2","T3","T4"), col = c(1,2,3,4), lty = 1)

#=================================================================================================
# 7 
# Zavislost Gamma na S

T1 <- 1
T2 <- 0.5
T3 <- 0.1
T4 <- 0.01

S <- seq(0,100, by = 1)

d1 <- rep(0,101)
d2 <- rep(0,101)
d3 <- rep(0,101)
d4 <- rep(0,101)

for (i in 1: 101){
  d1[i] <- Gama(S[i], K, T1, r, sigma)
  d2[i] <- Gama(S[i], K, T2, r, sigma)
  d3[i] <- Gama(S[i], K, T3, r, sigma)
  d4[i] <- Gama(S[i], K, T4, r, sigma)
  
}

plot(S, d1, type="l", xlab = "S", ylab = "Gama", xlim = c(0,100), ylim = c(0,0.3))
lines(S, d2, col = 2)
lines(S, d3, col = 3)
lines(S, d4, col = 4)
legend("topright", c("T1","T2","T3","T4"), col = c(1,2,3,4), lty = 1)

# ===========================================================================================
# 8
# Priebeh Delta call

S <- c(60, 40, 50)
t <- seq(0, 1, by = 0.01)
K <- 50
sigma <- 0.3
r <- 0

D.ITM <- rep(0, 101)
D.OTM <- rep(0, 101)
D.ATM <- rep(0, 101)

for (i in 1:101){
  D.ITM[i] <- Delta(S[1], K, t[i], r, sigma) 
  D.OTM[i] <- Delta(S[2], K, t[i], r, sigma) 
  D.ATM[i] <- Delta(S[3], K, t[i], r, sigma) 
}

plot(t, D.ITM, type = "l", ylim = c(0, 1.2), col = 1, ylab = "Delta")
lines(t, D.OTM, col = 2)
lines(t, D.ATM, col = 3)
legend("topright", c("ITM","OTM","ATM"), col = c(1,2,3), lty = 1, cex = 0.6 )

# ===========================================================================================
# 9
# Priebeh Gamma call

S <- c(60, 40, 50)
t <- seq(0, 1, by = 0.01)
K <- 50
sigma <- 0.3
r <- 0

G.ITM <- rep(0, 101)
G.OTM <- rep(0, 101)
G.ATM <- rep(0, 101)

for (i in 1:101){
  G.ITM[i] <- Gama(S[1], K, t[i], r, sigma) 
  G.OTM[i] <- Gama(S[2], K, t[i], r, sigma) 
  G.ATM[i] <- Gama(S[3], K, t[i], r, sigma) 
}

plot(t, G.ITM, type = "l", ylim = c(0, 0.3), col = 1, ylab = "Gamma")
lines(t, G.OTM, col = 2)
lines(t, G.ATM, col = 3)
legend("topright", c("ITM","OTM","ATM"), col = c(1,2,3), lty = 1 )

# ============================================================================================
# 10
# Cas poslednej navstevy pociatku

dt = 0.001
t = seq(0, 1, by = dt)
W0 = 0

W <- generuj.Wp (t = t, dt = dt, W0 = W0)
plot(t, W, type="l")
abline(h = W0, lty = 2)

posl.navsteva <- function(W){
  max.T <- NA
  for(i in 1:1000){
    if ((W[i] < W0 & W[i+1] > W0) | (W[i] > W0 & W[i+1] < W0)){
      max.T <- t[i]
    }
    if (is.na(max.T)) {max.T <- W0}
  }
  return(max.T)
}

max.T <- posl.navsteva(W)
abline(v = max.T, col = "red")

# ===========================================================================================
# 11

posl.navsteva2 <- function(){
  generuj.Wp <- function(mu = 0, sigma = 1, dt = 0.001, W0 = 0, t = seq(0, 1, by = dt)) {
    N <- length(t)
    dW <- rnorm(N - 1, mu, sigma) * sqrt(dt)
    W <- cumsum(c(W0, dW))
    return(W)
  }
  
  W <- generuj.Wp (t = t, dt = dt, W0 = W0)
  max.T <- NA
  for(i in 1:1000){
    if ((W[i] < W0 & W[i+1] > W0) | (W[i] > W0 & W[i+1] < W0)){
      max.T <- t[i]
    }
    if (is.na(max.T)) {max.T <- W0}
  }
  return(max.T)
}

n <- c(10, 100, 1000)

for (i in 1:3){
  maxT <- rep(0,n[i])
  for (j in 1:n[i]){
    maxT[j] <- posl.navsteva2()
  }
  hist(maxT,breaks = 100, main = c("Histogram pre n =", n[i]), plot = T)
  
  x <- seq(0,1,length.out = n[i])
  zakon <- n[i]/100 * 1/(pi*sqrt(x*(1-x)))
  lines(x, zakon, type = "l", col = "red")
}

# ===================================================================================
#12 Cash-or-Nothing

X0 <- 100 # cena akcie
K <- 50   # realizacni cena
Q <- 110  # vyplata

sigma <- 0.2
r <- 0
N <- 10000  # pocet trajektorii gBp

generuj.gBp <- function(X0, r, sigma, dt = 0.001, t = seq(0, 1, by = dt)){
  n <- length (t)
  dW <- rnorm (n - 1) * sqrt (dt)
  dX <- 1 + r * dt + sigma * dW
  X <- cumprod (c (X0 , dX))
}

cashornothing <- function(Q, X0, K, r, sigma, N, dt = 0.001, t= seq(0, 1, by = dt)){
  con <- rep(0,N)
  for (k in 1:N) {     
    W <- generuj.gBp (X0=X0,-(sigma^2)/2, sigma=sigma, dt = dt, t = t)
    if (W[length(W)] > K){con[k] <- 1}
  }
  Pst = sum(con)/N
  Ocekavana = Pst*Q
  tt <- t[length(t)]
  d2 <- (log (X0/K) + (r - sigma^2 / 2) * T) / sigma / sqrt (tt)
  Teoreticka <- exp(0) * Q * pnorm(d2)
  return(c(Ocekavana,Teoreticka))
}

cashornothing(Q, X0, K, r, sigma, N)

#==================================================================================
#13 
# Vasickuv model
# OUM - okamzita urokova mira

a <- 0.1
b <- 0.1
N <- 10   # pocet trajektorii
r0 <- 0.05
sigma <- 0.3
t <- seq(0, 1, by = dt)
dt <- 0.001

Vasicek <- function(N, a, b, t, dt, r0, sigma){
    OUM <- function(a=a, b=b, t=t, dt=dt, r0=r0, sigma=sigma){
    n <- length(t)
    r <- rep(0, n)
    r[1] <- r0
    for (k in 1:(n - 1)) {
      dW <- rnorm(1) * sqrt(dt)
      dr <- a * (b - r[k]) * dt + sigma * dW
      r[k + 1] <- r[k] + dr
    }
    return(r)
    }
    
    M <- sapply (1:N , function (k){
      OUM(a=a, b=b, t=t, dt=dt, r0=r0, sigma=sigma)
    })
    matplot(t, M, type = "l", lty = 1, col = "#ccaaee", xlab = "t", ylab = "r", 
            main = "Trajektoria r")
}

Vasicek(N, a, b, t, dt, r0, sigma)

#==================================================================================
#14 
# CIR model

a <- 0.1
b <- 0.1
N <- 10   # pocet trajektorii
r0 <- 0.05
sigma <- 0.3
t <- seq(0, 1, by = dt)
dt <- 0.001

CIR <- function(N, a, b, t, dt, r0, sigma){
  generuj.CIR <- function(a=a, b=b, t=t, dt=dt, r0=r0, sigma=sigma){
    n <- length(t)
    r <- rep(0, n)
    r[1] <- r0
    for (k in 1:(n - 1)) {
      dW <- rnorm(1) * sqrt(dt)
      if (r[k] < 0){r[k] <- 0} 
      dr <- a * (b - r[k]) * dt + sigma * sqrt(r[k]) * dW
      r[k + 1] <- r[k] + dr
    }
    return(r)
   }
  
  M <- sapply (1:N , function (k){
    generuj.CIR(a=a, b=b, t=t, dt=dt, r0=r0, sigma=sigma)
  })
  matplot(t, M, type = "l", lty = 1, col = "#ccaaee", xlab = "t", ylab = "r", 
          main = "Trajektorie r")
}

CIR(N, a, b, t, dt, r0, sigma)

# ====================================================================================
#15
# Implikovana volatilita

library(readxl)
FB <- read_excel("FB.xlsx")

sigmaimpl <- rep(0,length(FB$K))

for(i in 1:170){
  fn <- function (u) {
    call.opce (FB$K[i], FB$S[1], u, FB$r[i]/100, FB$t[i]) - FB$V[i]
  }
  reseni<- uniroot (fn,c(0,2))
  sigmaimpl[i]<-reseni$root
}

plot(FB$K[1:25],sigmaimpl[1:25],type="l")       # t = 0,013
plot(FB$K[26:49],sigmaimpl[26:49],type="l")     # t = 0,11
plot(FB$K[50:74],sigmaimpl[50:74],type="l")     # t = 0,186
plot(FB$K[75:95],sigmaimpl[75:95],type="l")     # t = 0,263
plot(FB$K[96:120],sigmaimpl[96:120],type="l")   # t = 0,435
plot(FB$K[121:145],sigmaimpl[121:145],type="l") # t = 1,01
plot(FB$K[146:170],sigmaimpl[146:170],type="l") # t = 2

# Plocha implikovanej volatility
plot3d(FB$K,sigmaimpl,FB$t) 

#======================================================================================
# 16
# Aritmeticky priemer WP

dt <- 0.001
t <- seq(0, 1, by = dt)
X0 <- 1
mu <- 0.5
sigma <- 0.2

X <- generuj.Wp()
plot(t, X, type = "l", lty = 1, main = "trajektorie WP", xlab = "t",
     ylab = "X")
abline(h = X0 , lty = 2)

AP <- function(t, dt, X0) {
  n <- length(t)
  Pr <- rep(0,n)
  for (i in 1:n) {
    Pr[i] <- (1/i) * sum(X[1:i])
  }
  return(Pr)
}

priemer <- AP(t, dt, X0)
plot(t, X, type = "l", xlab = "t", ylab = "X", main="Trajektoria aritmetickeho priemeru")
lines(t, priemer, type = "l", col = "red")
abline(h = X0, lty = 2)

#======================================================================================
# 17
# Azijska opcia

X0<- 100
r <- 0.05
sigma <- 0.3
t <- 1

generuj.gBpPrumer <- function (X0, r, sigma, dt = 0.001, t = seq(0, 1, by=dt)){
  n <- length (t)
  dW <- rnorm (n - 1) * sqrt (dt)
  dX <- 1+r * dt + sigma * dW
  X <- cumprod(c (X0 , dX))
  Z=cumsum(X)/(1:n)
  if(X[n]-Z[n]>0){C=X[n]-Z[n]}
  else{C=0}
  return(C)
}

asian.option <-function(X0,sigma, N = 1000, dt = 0.001, t = seq(0, 1, by=dt)){
  suma<-rep(0,N)
  for (k in 1:N) {
    suma[k]<-generuj.gBpPrumer(X0 = X0, -(sigma^2)/2, sigma = sigma, dt = dt, t = t)
  }
  return(sum(suma)/N)
}

asian.option(X0,sigma)