# ------------------------------------------------------
# Implémentation de la méthode de Newton-Raphson pour 
# le calcul de la volatilité implicite selon Black-Scholes
# ------------------------------------------------------

# Fonction de prix d'option selon Black-Scholes avec dividendes
black_scholes_call <- function(S, K, T, r, q, sigma) {
  # Calcul des paramètres d1 et d2
  d1 <- (log(S/K) + (r - q + sigma^2/2)*T) / (sigma * sqrt(T))
  d2 <- d1 - sigma * sqrt(T)
  
  # Prix du call selon la formule avec dividendes
  return(S * exp(-q*T) * pnorm(d1) - K * exp(-r*T) * pnorm(d2))
}

# Vega de l'option (dérivée du prix par rapport à sigma)
vega <- function(S, K, T, r, q, sigma) {
  d1 <- (log(S/K) + (r - q + sigma^2/2)*T) / (sigma * sqrt(T))
  # Formule du vega: S * e^(-q*T) * n(d1) * sqrt(T)
  # où n(d1) est la fonction de densité de la loi normale standard
  return(S * exp(-q*T) * dnorm(d1) * sqrt(T))
}


# Méthode de Newton pour trouver la volatilité implicite
implied_volatility_newton <- function(market_price, S, K, T, r, q, 
                                      initial_vol = 0.2,
                                      max_iterations = 100, 
                                      precision = 1e-6) {
  # Initialisation de la volatilité
  sigma <- initial_vol
  
  # Pour tracer la convergence
  iterations <- c(0)
  sigmas <- c(sigma)
  prices <- c(black_scholes_call(S, K, T, r, q, sigma))
  diffs <- c(prices[1] - market_price)
  vegas <- c(vega(S, K, T, r, q, sigma))
  adjustments <- c(NA)  # Premier ajustement N/A
  
  # Boucle principale de Newton-Raphson
  for (i in 1:max_iterations) {
    # Calcul du prix avec la volatilité courante
    price <- black_scholes_call(S, K, T, r, q, sigma)
    price_diff <- price - market_price
    
    # Vérification de la convergence
    if (abs(price_diff) < precision) {
      cat(sprintf("Convergence en %d itérations\n", i))
      break
    }
    
    # Calcul du vega (dérivée du prix par rapport à sigma)
    v <- vega(S, K, T, r, q, sigma)
    
    # Éviter la division par zéro
    if (v == 0) {
      sigma <- sigma * 1.5
      next
    }
    
    # Méthode de Newton: x_n+1 = x_n - f(x_n)/f'(x_n)
    adjustment <- price_diff / v
    sigma_new <- sigma - adjustment
    
    # Borner la volatilité pour éviter les valeurs négatives ou extrêmes
    sigma_new <- max(0.001, min(sigma_new, 5))
    
    # Enregistrement pour le suivi
    iterations <- c(iterations, i)
    sigmas <- c(sigmas, sigma_new)
    prices <- c(prices, price)
    diffs <- c(diffs, price_diff)
    vegas <- c(vegas, v)
    adjustments <- c(adjustments, adjustment)
    
    # Mise à jour pour la prochaine itération
    sigma <- sigma_new
  }
  
  # Retourner les résultats
  return(list(implied_vol = sigma,
              iterations = iterations,
              sigmas = sigmas,
              prices = prices,
              diffs = diffs,
              vegas = vegas,
              adjustments = adjustments))
}

# ------------------------------------------------------
# Exemple d'application sur le cas du CAC40
# ------------------------------------------------------

# Paramètres de l'option comme présenté dans le PDF
S <- 7900           # Prix du sous-jacent (CAC40)
K <- 8100           # Prix d'exercice
T <- 0.25           # Temps à maturité (3 mois)
r <- 0.03           # Taux sans risque (3%)
q <- 0.02           # Taux de rendement du dividende (2%)
true_vol <- 0.19    # Volatilité historique
market_price <- 200 # Prix de marché du call observé

# Calculer le prix avec la volatilité historique
initial_price <- black_scholes_call(S, K, T, r, q, true_vol)
cat(sprintf("Prix avec volatilité historique (%0.2f%%): %0.2f???\n", 
            true_vol*100, initial_price))

# Trouver la volatilité implicite avec différentes valeurs initiales
initial_vols <- c(0.1, 0.2, 0.3)
results <- list()

# Calcul pour chaque valeur initiale
for (i in 1:length(initial_vols)) {
  initial_vol <- initial_vols[i]
  result <- implied_volatility_newton(market_price, S, K, T, r, q, initial_vol)
  results[[i]] <- result
  
  cat(sprintf("Vol. initiale: %0.2f%%, Vol. implicite trouvée: %0.4f%%, Itérations: %d\n", 
              initial_vol*100, result$implied_vol*100, length(result$iterations)-1))
}

# ------------------------------------------------------
# Visualisation des résultats
# ------------------------------------------------------

# Configurer un graphique à deux panneaux
par(mfrow = c(2, 1), mar = c(4, 4, 3, 1))
colors <- c("blue", "green", "purple")

# Premier graphique pour la convergence
plot(NA, xlim = c(0, 10), ylim = c(0.05, 0.35),
     xlab = "Itérations", ylab = "Volatilité estimée",
     main = "Convergence de la méthode de Newton pour différentes valeurs initiales")
abline(h = results[[1]]$implied_vol, col = "red", lty = 2, lwd = 2)
legend("topright", legend = c(paste("Vol. initiale =", initial_vols*100, "%"), "Vol. implicite"),
       col = c(colors, "red"), lty = c(rep(1, 3), 2), lwd = 2)

# Tracer la convergence pour chaque valeur initiale
for (i in 1:length(results)) {
  result <- results[[i]]
  lines(result$iterations, result$sigmas, type = "o", col = colors[i], lwd = 2)
}

# Deuxième graphique pour la fonction objectif (prix_BS - prix_marché)
sigmas_range <- seq(0.05, 0.4, length.out = 100)
price_diffs <- sapply(sigmas_range, function(sigma) {
  black_scholes_call(S, K, T, r, q, sigma) - market_price
})

plot(sigmas_range, price_diffs, type = "l", col = "blue", lwd = 2,
     xlab = "Volatilité", ylab = "Prix BS - Prix marché",
     main = "Recherche du zéro de la fonction f(??) = Prix_BS(??) - Prix_marché")
abline(h = 0, col = "red", lty = 2, lwd = 2)
abline(v = results[[1]]$implied_vol, col = "green", lty = 2, lwd = 2)
legend("topleft", legend = c("f(??) = Prix_BS(??) - Prix_marché", "Zéro de la fonction", "Volatilité implicite"),
       col = c("blue", "red", "green"), lty = c(1, 2, 2), lwd = 2)

# Création d'un tableau des itérations
# Sélectionner le résultat avec vol. initiale = 0.2 pour le tableau
result <- results[[2]]  # Correspondant à initial_vol = 0.2
n_iter <- length(result$iterations) - 1  # Nombre d'itérations réelles

# Créer un data frame pour faciliter l'affichage
iterations_df <- data.frame(
  Iteration = 0:n_iter,
  Volatilite = round(result$sigmas * 100, 2),
  Prix = round(result$prices, 2),
  Difference = round(result$diffs, 2),
  Vega = round(result$vegas, 2),
  Ajustement = c(NA, round(result$adjustments[-1] * 100, 4))  # En points de pourcentage
)

# Afficher le tableau
cat("\nTableau des itérations de Newton-Raphson:\n")
print(iterations_df)

# Afficher les résultats finaux
cat("\nRésultats finaux:\n")
cat(sprintf("Volatilité historique: %.2f%%\n", true_vol * 100))
cat(sprintf("Volatilité implicite: %.2f%%\n", result$implied_vol * 100))
cat(sprintf("Décote de volatilité: %.2f points\n", (result$implied_vol - true_vol) * 100))
cat(sprintf("Décote relative: %.2f%%\n", (result$implied_vol/true_vol - 1) * 100))