Libraries

# install.packages("dplyr")
# install.packages("quantmod")
# install.packages("moments")
# install.packages("ggplot2")
# install.packages("gridExtra")
# install.packages("fitdistrplus")
# install.packages("marida")
# install.packages("nortest")
# install.packages("copula")
# install.packages("DT")
# install.packages("MVN")

suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(quantmod))
suppressPackageStartupMessages(library(moments))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(gridExtra))
suppressPackageStartupMessages(library(fitdistrplus))
suppressPackageStartupMessages(library(nortest))
suppressPackageStartupMessages(library(copula))
suppressPackageStartupMessages(library(plotly))
suppressPackageStartupMessages(library(DT))
suppressPackageStartupMessages(library(MVN))

Data Collecting

nvda <- getSymbols("NVDA", src = "yahoo", from = "2020-10-30", to = "2025-10-30", auto.assign = FALSE)[,6] # "yyyy-mm-dd"
orcl <- getSymbols("ORCL", src = "yahoo", from = "2020-10-30", to = "2025-10-30", auto.assign = FALSE)[,6] # "yyyy-mm-dd"

Data Plotting

# Convert xts objects to data frames for plotly
nvda_df <- data.frame(Date = index(nvda), Price = as.numeric(nvda))
colnames(nvda_df)[2] <- "NVDA_Price"

orcl_df <- data.frame(Date = index(orcl), Price = as.numeric(orcl))
colnames(orcl_df)[2] <- "ORCL_Price"

# Plot 1: NVDA 
nvda_plot <- plot_ly(nvda_df, x = ~Date, y = ~NVDA_Price,
                     type = 'scatter', mode = 'lines',
                     line = list(color = 'darkgreen', width = 2),
                     name = 'NVDA',
                     hovertemplate = 'Date: %{x}<br>Price: $%{y:.2f}<extra></extra>') %>%
  layout(title = list(text = "NVDA Stock Price (2020-2025)", x = 0.5),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Price (USD)"))

# Plot 2: ORCL 
orcl_plot <- plot_ly(orcl_df, x = ~Date, y = ~ORCL_Price,
                     type = 'scatter', mode = 'lines',
                     line = list(color = 'darkblue', width = 2),
                     name = 'ORCL',
                     hovertemplate = 'Date: %{x}<br>Price: $%{y:.2f}<extra></extra>') %>%
  layout(title = list(text = "ORCL Stock Price (2020-2025)", x = 0.5),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Price (USD)"))

# Plot 3: combined NVDA & ORCL 
combined_plot <- plot_ly() %>%
  add_trace(data = nvda_df, x = ~Date, y = ~NVDA_Price,
            type = 'scatter', mode = 'lines',
            line = list(color = 'darkgreen', width = 2),
            name = 'NVDA',
            hovertemplate = 'Date: %{x}<br>NVDA: $%{y:.2f}<extra></extra>') %>%
  add_trace(data = orcl_df, x = ~Date, y = ~ORCL_Price,
            type = 'scatter', mode = 'lines',
            line = list(color = 'darkblue', width = 2),
            name = 'ORCL',
            hovertemplate = 'Date: %{x}<br>ORCL: $%{y:.2f}<extra></extra>',
            yaxis = "y") %>%
  layout(title = list(text = "NVDA vs ORCL Stock Price (2020-2025)", x = 0.5),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Price (USD)",
                      rangemode = "tozero"),
         legend = list(x = 0.02, y = 0.98,
                       bgcolor = 'rgba(255,255,255,0.8)'),
         hovermode = "x unified")

# Display plots
nvda_plot
orcl_plot
combined_plot

Obliczanie stop strat

nvda_loss <- round(c(NA, -diff(as.numeric(Cl(to.monthly(nvda)))) / head(as.numeric(Cl(to.monthly(nvda))), -1) * 100), 2)
orcl_loss <- round(c(NA, -diff(as.numeric(Cl(to.monthly(orcl)))) / head(as.numeric(Cl(to.monthly(orcl))), -1) * 100), 2)
nvda_monthly <- as.numeric(Cl(to.monthly(nvda)))
orcl_monthly <- as.numeric(Cl(to.monthly(orcl)))

loss_data <- data.frame(
  Data = index(to.monthly(nvda)),
  nvda_close = nvda_monthly,
  nvda_loss = nvda_loss,
  orcl_close = orcl_monthly,
  orcl_loss = orcl_loss
)
loss_data <- loss_data %>% filter(!is.na(nvda_loss))

Statystyki opisowe

summary_stats <- data.frame(
  Statistic = c("N", "Średnia", "Mediana", "Wariancja", "Odchylenie standardowe", "Skośność", "Kurtoza"),
  NVDA = c(
    length(na.omit(nvda_loss)),
    mean(nvda_loss, na.rm = TRUE),
    median(nvda_loss, na.rm = TRUE),
    var(nvda_loss, na.rm = TRUE),
    sd(nvda_loss, na.rm = TRUE),
    skewness(nvda_loss, na.rm = TRUE),
    kurtosis(nvda_loss, na.rm = TRUE) - 3
  ),
  ORCL = c(
    length(na.omit(orcl_loss)),
    mean(orcl_loss, na.rm = TRUE),
    median(orcl_loss, na.rm = TRUE),
    var(orcl_loss, na.rm = TRUE),
    sd(orcl_loss, na.rm = TRUE),
    skewness(orcl_loss, na.rm = TRUE),
    kurtosis(orcl_loss, na.rm = TRUE) - 3
  )
)

summary_stats
##                Statistic        NVDA        ORCL
## 1                      N  60.0000000  60.0000000
## 2                Średnia  -5.8018333  -3.3126667
## 3                Mediana  -5.7550000  -2.2900000
## 4              Wariancja 210.5815644 108.8478402
## 5 Odchylenie standardowe  14.5114288  10.4330168
## 6               Skośność   0.1207483  -0.4984035
## 7                Kurtoza  -0.3362411   0.1624986

Testy normalności

shapiro.test(loss_data$nvda_loss)
## 
##  Shapiro-Wilk normality test
## 
## data:  loss_data$nvda_loss
## W = 0.99196, p-value = 0.9627
ad.test(loss_data$nvda_loss)
## 
##  Anderson-Darling normality test
## 
## data:  loss_data$nvda_loss
## A = 0.13879, p-value = 0.9737
shapiro.test(loss_data$orcl_loss)
## 
##  Shapiro-Wilk normality test
## 
## data:  loss_data$orcl_loss
## W = 0.97655, p-value = 0.3006
ad.test(loss_data$orcl_loss)
## 
##  Anderson-Darling normality test
## 
## data:  loss_data$orcl_loss
## A = 0.46017, p-value = 0.2523

Histogramy

par(mfrow = c(1, 2), mar = c(5, 4, 3, 1))

# Histogram NVDA
h1 <- hist(loss_data$nvda_loss, 
     breaks = 10,
     main = "Rozkład L_nvda",
     xlab = "L_nvda",
     ylab = "Gęstość",
     col = rgb(0.5, 0.7, 0.9, 0.6),
     border = "white",
     freq = FALSE,
     las = 1,
     xaxt = "n",
     cex.main = 1.3,
     cex.lab = 1.1)
axis(1, at = h1$breaks, labels = round(h1$breaks, 1), las = 2, cex.axis = 0.8)
curve(dnorm(x, mean = mean(loss_data$nvda_loss), sd = sd(loss_data$nvda_loss)),
      add = TRUE, col = "darkblue", lwd = 2.5)

# Histogram ORCL
h2 <- hist(loss_data$orcl_loss,
     breaks = 10,
     main = "Rozkład L_orcl",
     xlab = "L_orcl",
     ylab = "Gęstość",
     col = rgb(0.5, 0.7, 0.9, 0.6),
     border = "white",
     freq = FALSE,
     las = 1,
     xaxt = "n",
     cex.main = 1.3,
     cex.lab = 1.1)
axis(1, at = h2$breaks, labels = round(h2$breaks, 1), las = 2, cex.axis = 0.8)
curve(dnorm(x, mean = mean(loss_data$orcl_loss), sd = sd(loss_data$orcl_loss)),
      add = TRUE, col = "darkblue", lwd = 2.5)

par(mfrow = c(1, 1))

Parametry rozkładu normalnego

mean_nvda <- mean(loss_data$nvda_loss)
sd_nvda <- sd(loss_data$nvda_loss)

mean_orcl <- mean(loss_data$orcl_loss)
sd_orcl <- sd(loss_data$orcl_loss)

print(paste("Średnia NVDA:", round(mean_nvda, 2), "%"))
## [1] "Średnia NVDA: -5.8 %"
print(paste("Odchylenie standardowe NVDA:", round(sd_nvda, 2), "%"))
## [1] "Odchylenie standardowe NVDA: 14.51 %"
print(paste("Średnia ORCL:", round(mean_orcl, 2), "%"))
## [1] "Średnia ORCL: -3.31 %"
print(paste("Odchylenie standardowe ORCL:", round(sd_orcl, 2), "%"))
## [1] "Odchylenie standardowe ORCL: 10.43 %"

VaR

empirical_nvda_VaR <- quantile(loss_data$nvda_loss, probs = 0.95)
empirical_orcl_VaR <- quantile(loss_data$orcl_loss, probs = 0.95)

theoretical_nvda_VaR <- mean_nvda + sd_nvda * qnorm(0.95,mean=0,sd=1)
theoretical_orcl_VaR <- mean_orcl + sd_orcl * qnorm(0.95,mean=0,sd=1)

print(paste("Empiryczny VaR NVDA (95%):", round(empirical_nvda_VaR, 2), "%"))
## [1] "Empiryczny VaR NVDA (95%): 16.99 %"
print(paste("Teoretyczny VaR NVDA (95%):", round(theoretical_nvda_VaR, 2), "%"))
## [1] "Teoretyczny VaR NVDA (95%): 18.07 %"
print(paste("Empiryczny VaR ORCL (95%):", round(empirical_orcl_VaR, 2), "%"))
## [1] "Empiryczny VaR ORCL (95%): 10.98 %"
print(paste("Teoretyczny VaR ORCL (95%):", round(theoretical_orcl_VaR, 2), "%"))
## [1] "Teoretyczny VaR ORCL (95%): 13.85 %"

Expected Shortfall

ES_empirical_nvda <- mean(loss_data$nvda_loss[loss_data$nvda_loss >= empirical_nvda_VaR])
ES_empirical_orcl <- mean(loss_data$orcl_loss[loss_data$orcl_loss >= empirical_orcl_VaR])
         
ES_theoretical_nvda1 <- mean_nvda + sd_nvda * (dnorm(qnorm(0.95)) / (1 - 0.95))
ES_theoretical_orcl1 <- mean_orcl + sd_orcl * (dnorm(qnorm(0.95)) / (1 - 0.95))
                                            
print(paste("NVIDIA - Wartość wyznaczona za pomocą MPWL:", round(ES_empirical_nvda, 2), "%"))
## [1] "NVIDIA - Wartość wyznaczona za pomocą MPWL: 23.46 %"
print(paste("Oracle - Wartość wyznaczona za pomocą MPWL:", round(ES_empirical_orcl, 2), "%"))
## [1] "Oracle - Wartość wyznaczona za pomocą MPWL: 15.16 %"
print(paste("NVIDIA - Wartość wyznaczona za pomocą wzoru analitycznego:", round(ES_theoretical_nvda1, 2), "%"))
## [1] "NVIDIA - Wartość wyznaczona za pomocą wzoru analitycznego: 24.13 %"
print(paste("Oracle - Wartość wyznaczona za pomocą wzoru analitycznego:", round(ES_theoretical_orcl1, 2), "%"))
## [1] "Oracle - Wartość wyznaczona za pomocą wzoru analitycznego: 18.21 %"

VaR & ES Vizualizacja

d_orcl <- density(loss_data$orcl_loss)
d_nvda <- density(loss_data$nvda_loss)

# ORCL - zaawansowane śledzenie kursora
orcl_plot <- plot_ly() %>%
  # Gęstość empiryczna
  add_trace(x = d_orcl$x, y = d_orcl$y,
            type = 'scatter', mode = 'lines',
            line = list(color = 'darkblue', width = 2),
            fill = 'tozeroy',
            fillcolor = 'rgba(0,0,255,0.2)',
            name = 'Gęstość empiryczna') %>%
  
  # Normalna krzywa
  add_trace(x = seq(min(d_orcl$x), max(d_orcl$x), length=200),
            y = dnorm(seq(min(d_orcl$x), max(d_orcl$x), length=200), 
                     mean=-3.31, sd=10.43),
            type = 'scatter', mode = 'lines',
            line = list(color = 'black', width=2),
            name = 'Rozkład normalny') %>%
  
  # Empiryczny ES (punkt)
  add_trace(x = ES_empirical_orcl, y = 0,
            type = 'scatter', mode = 'markers',
            marker = list(color='red', size=15),
            name = 'ES Empiryczny') %>%
  
  # Teoretyczny ES (punkt)
  add_trace(x = ES_theoretical_orcl1, y = 0,
            type = 'scatter', mode = 'markers',
            marker = list(color='orange', size=15),
            name = 'ES Teoretyczny') %>%
  
  layout(
    title = list(text = "VaR & ES (α=0.95): ORCL", x = 0.5),
    xaxis = list(
      title = "Strata (%)",
      zeroline = FALSE,
      showspikes = TRUE,
      spikemode = "across+toaxis",
      spikesnap = "cursor",
      spikedash = "dash",
      spikecolor = "#888",
      spikethickness = 1,
      showline = TRUE,
      linecolor = "black"
    ),
    yaxis = list(
      title = "Gęstość",
      zeroline = FALSE,
      showspikes = TRUE,
      spikemode = "across+toaxis",
      spikesnap = "cursor",
      spikedash = "dash",
      spikecolor = "#888",
      spikethickness = 1,
      showline = TRUE,
      linecolor = "black"
    ),
    hovermode = "x",
    showlegend = TRUE,
    dragmode = "zoom",
    plot_bgcolor = 'white'
  ) %>%
  config(
    displayModeBar = TRUE,
    scrollZoom = TRUE,
    modeBarButtonsToAdd = list("drawline", "drawopenpath", "eraseshape")
  )

# NVDA - zaawansowane śledzenie kursora
nvda_plot <- plot_ly() %>%
  # Gęstość empiryczna
  add_trace(x = d_nvda$x, y = d_nvda$y,
            type = 'scatter', mode = 'lines',
            line = list(color = 'darkgreen', width = 2),
            fill = 'tozeroy',
            fillcolor = 'rgba(0,255,0,0.2)',
            name = 'Gęstość empiryczna') %>%
  
  # Normalna krzywa
  add_trace(x = seq(min(d_nvda$x), max(d_nvda$x), length=200),
            y = dnorm(seq(min(d_nvda$x), max(d_nvda$x), length=200), 
                     mean=-5.8, sd=14.51),
            type = 'scatter', mode = 'lines',
            line = list(color = 'black', width=2),
            name = 'Rozkład normalny') %>%
  
  # Empiryczny ES (punkt)
  add_trace(x = ES_empirical_nvda, y = 0,
            type = 'scatter', mode = 'markers',
            marker = list(color='red', size=15),
            name = 'ES Empiryczny') %>%
  
  # Teoretyczny ES (punkt)
  add_trace(x = ES_theoretical_nvda1, y = 0,
            type = 'scatter', mode = 'markers',
            marker = list(color='orange', size=15),
            name = 'ES Teoretyczny') %>%
  
  layout(
    title = list(text = "VaR & ES (α=0.95): NVDA", x = 0.5),
    xaxis = list(
      title = "Strata (%)",
      zeroline = FALSE,
      showspikes = TRUE,
      spikemode = "across+toaxis",
      spikesnap = "cursor",
      spikedash = "dash",
      spikecolor = "#888",
      spikethickness = 1,
      showline = TRUE,
      linecolor = "black"
    ),
    yaxis = list(
      title = "Gęstość", 
      range = c(0, 0.03),
      zeroline = FALSE,
      showspikes = TRUE,
      spikemode = "across+toaxis",
      spikesnap = "cursor",
      spikedash = "dash",
      spikecolor = "#888",
      spikethickness = 1,
      showline = TRUE,
      linecolor = "black"
    ),
    hovermode = "x",
    showlegend = TRUE,
    dragmode = "zoom",
    plot_bgcolor = 'white'
  ) %>%
  config(
    displayModeBar = TRUE,
    scrollZoom = TRUE,
    modeBarButtonsToAdd = list("drawline", "drawopenpath", "eraseshape")
  )

# Wyświetl
orcl_plot
nvda_plot

Dane dla portfela dwóch akcji

laczny <- data.frame(nvda = loss_data$nvda_loss, orcl = loss_data$orcl_loss)

Mardina test

rezultat_mardia <- mvn(laczny, mvn_test = "mardia")
rezultat_mardia
## $multivariate_normality
##              Test Statistic p.value     Method          MVN
## 1 Mardia Skewness    10.539   0.032 asymptotic ✗ Not normal
## 2 Mardia Kurtosis    -0.452   0.652 asymptotic     ✓ Normal
## 
## $univariate_normality
##               Test Variable Statistic p.value Normality
## 1 Anderson-Darling     nvda     0.139   0.974  ✓ Normal
## 2 Anderson-Darling     orcl     0.460   0.252  ✓ Normal
## 
## $descriptives
##   Variable  n   Mean Std.Dev Median    Min   Max    25th  75th   Skew Kurtosis
## 1     nvda 60 -5.802  14.511 -5.755 -36.34 32.03 -14.722 2.695  0.121    2.664
## 2     orcl 60 -3.313  10.433 -2.290 -32.08 17.64 -10.220 2.475 -0.498    3.162
## 
## $data
##      nvda   orcl
## 1   -6.92  -2.87
## 2    2.56 -12.08
## 3    0.50   6.23
## 4   -5.58  -6.75
## 5    2.64  -8.77
## 6  -12.45  -8.48
## 7   -8.23  -3.89
## 8  -23.16   1.14
## 9    2.52 -12.36
## 10 -14.82  -2.28
## 11   7.46   2.24
## 12 -23.42 -10.49
## 13 -27.81   5.42
## 14   9.98   3.89
## 15  16.75   6.59
## 16   0.41   6.39
## 17 -11.92  -8.90
## 18  32.03  10.93
## 19  -0.67   2.02
## 20  18.80   2.85
## 21 -19.82 -11.91
## 22  16.90   4.74
## 23  19.55  17.64
## 24 -11.19 -28.49
## 25 -25.42  -6.35
## 26  13.64   1.55
## 27 -33.69  -8.63
## 28 -18.83   1.20
## 29 -19.67  -6.32
## 30   0.10  -2.36
## 31 -36.34 -11.85
## 32 -11.82 -12.41
## 33 -10.47   1.22
## 34  -5.62  -2.70
## 35  11.86  12.02
## 36   6.25   2.02
## 37 -14.69 -12.39
## 38  -5.89   9.28
## 39 -24.24  -6.36
## 40 -28.58   0.02
## 41 -14.22 -12.47
## 42   4.38   9.15
## 43 -26.89  -3.02
## 44 -12.69 -20.49
## 45   5.28   0.96
## 46  -2.01  -1.32
## 47  -1.74 -20.60
## 48  -9.32   1.28
## 49  -4.14 -10.13
## 50   2.86   9.85
## 51  10.59  -2.30
## 52  -4.04   2.35
## 53  13.23  15.81
## 54  -0.50  -1.01
## 55 -24.06 -17.63
## 56 -16.93 -32.08
## 57 -12.58 -16.32
## 58   2.07  10.89
## 59  -7.13 -24.37
## 60 -10.97   1.94
## 
## $subset
## NULL
## 
## $outlierMethod
## [1] "none"
## 
## attr(,"class")
## [1] "mvn"

Dopasowanie kopuly

# Przekształć dane na rozkłady jednostajne [0,1]
u_data <- pobs(as.matrix(laczny))

# A. Kopuły Archimedejskie
fit_clayton <- fitCopula(claytonCopula(), u_data, method='ml')
fit_frank <- fitCopula(frankCopula(), u_data, method='ml')
fit_gumbel <- fitCopula(gumbelCopula(), u_data, method='ml')

# B. Kopuła eliptyczna
fit_t <- fitCopula(tCopula(), u_data, method='ml')
fit_normal <- fitCopula(normalCopula(), u_data, method='ml')

# Stwórz tabelę porównawczą
comparison <- data.frame(
  Kopula = c("Clayton", "Frank", "Gumbel", "t-Student", "Normal"),
  LogLik = c(logLik(fit_clayton), logLik(fit_frank), 
             logLik(fit_gumbel), logLik(fit_t), logLik(fit_normal)),
  AIC = c(AIC(fit_clayton), AIC(fit_frank), 
          AIC(fit_gumbel), AIC(fit_t), AIC(fit_normal)),
  BIC = c(BIC(fit_clayton), BIC(fit_frank), 
          BIC(fit_gumbel), BIC(fit_t), BIC(fit_normal))
)

# Posortuj wyniki
comparison[order(comparison$AIC),]  # Najniższe AIC = najlepsze
##      Kopula    LogLik        AIC        BIC
## 3    Gumbel 12.682784 -23.365567 -21.271223
## 5    Normal 10.932616 -19.865232 -17.770887
## 2     Frank 10.883187 -19.766374 -17.672030
## 4 t-Student 10.930496 -17.860992 -13.672303
## 1   Clayton -1.205451   4.410902   6.505246
# Generate grid
grid_seq <- seq(0.01, 0.99, length.out = 30)
grid <- expand.grid(u = grid_seq, v = grid_seq)

# Calculate values
dens_vals <- dCopula(as.matrix(grid), fit_gumbel@copula)
dist_vals <- pCopula(as.matrix(grid), fit_gumbel@copula)

# Create matrices
dens_matrix <- matrix(dens_vals, 30, 30)
dist_matrix <- matrix(dist_vals, 30, 30)

# Distribution plot  
plot_ly(z = ~dist_matrix, type = "surface", colorscale = "Reds") %>%
  layout(scene = list(xaxis = list(title = "u"),
                      yaxis = list(title = "v"),
                      zaxis = list(title = "C(u,v)")))

Wizualizacja kopuli 2D & 3D

# 2D CONTOUR PLOT OF COPULA DENSITY

# 1. Grid
nx <- 100
u <- seq(0, 1, length.out = nx)
v <- seq(0, 1, length.out = nx)
grid <- expand.grid(u = u, v = v)

# 2. Copula density matrix
z <- matrix(
  dCopula(cbind(grid$u, grid$v), fit_gumbel@copula),
  nrow = nx, 
  ncol = nx
)

# 3. Plotly contour + heatmap
fig2D <- plot_ly(
  x = u, 
  y = v, 
  z = z,
  type = "contour",
  colorscale = "Viridis",
  contours = list(
    coloring = "heatmap",
    showlabels = TRUE
  )
) %>%
  layout(
    title = "2D Density Contour Plot – Gumbel Copula",
    xaxis = list(title = "u (NVDA)"),
    yaxis = list(title = "v (ORCL)")
  )

fig2D
# 3D
# 1. Define Gumbel copula
gumbel_theta <- 1.651
gumbel_cop <- gumbelCopula(param = gumbel_theta, dim = 2)

# 2. Grid
nx <- 100
u <- seq(0, 1, length.out = nx)
v <- seq(0, 1, length.out = nx)
grid <- expand.grid(u = u, v = v)

# 3. Copula density
z <- matrix(dCopula(cbind(grid$u, grid$v), gumbel_cop), nrow = nx, ncol = nx)

# 4. 3D Surface Plot
fig <- plot_ly(x = u, y = v, z = z) %>%
  add_surface(contours = list(
    z = list(show=TRUE, usecolormap=TRUE, highlightcolor="#ff0000", project=list(z=TRUE))
  )) %>%
  layout(
    title = "3D Density Surface – Gumbel Copula",
    scene = list(
      xaxis = list(title = "u (NVDA)"),
      yaxis = list(title = "v (ORCL)"),
      zaxis = list(title = "Density")
    )
  )

fig

Tabela & Wykres portfolio

# 1. Generowanie danych z kopuły (100 000 obserwacji)
cop_data <- rCopula(100000, fit_gumbel@copula)

# 2. Przekształcenie na rozkłady brzegowe
nvda_sim <- qnorm(cop_data[,1], mean_nvda, sd_nvda)  
orcl_sim <- qnorm(cop_data[,2], mean_orcl, sd_orcl)

# Wagi beta
beta_seq <- seq(0, 1, 0.1)
results <- data.frame(beta = beta_seq, VaR = NA, ES = NA)

for (i in 1:length(beta_seq)) {
  beta <- beta_seq[i]
  portfolio <- beta * nvda_sim + (1-beta) * orcl_sim
  
  VaR_95 <- quantile(portfolio, 0.95)
  ES_95 <- mean(portfolio[portfolio >= VaR_95])
  
  results$VaR[i] <- round(VaR_95, 6)
  results$ES[i] <- round(ES_95, 6)
}

# Tworzenie interaktywnych wykresów za pomocą Plotly

Sys.setlocale("LC_ALL", "Polish_Poland.UTF-8")
## [1] "LC_COLLATE=Polish_Poland.utf8;LC_CTYPE=Polish_Poland.utf8;LC_MONETARY=Polish_Poland.utf8;LC_NUMERIC=C;LC_TIME=Polish_Poland.utf8"
options(encoding = "UTF-8")

# Wykres 1: VaR (95%)
var_plot <- plot_ly(results, x = ~beta, y = ~VaR,
                    type = 'scatter', mode = 'lines+markers',
                    line = list(color = 'darkblue', width = 3),
                    marker = list(size = 8, color = 'darkblue'),
                    name = 'VaR',
                    hovertemplate = paste(
                      'Beta: %{x:.1f}<br>',
                      'VaR: %{y:.4f}<extra></extra>'
                    )) %>%
  layout(title = list(text = "Wartość Zagrożona (VaR 95%)", x = 0.5),
         xaxis = list(title = "Waga na NVDA (β)",
                      tickmode = 'array',
                      tickvals = beta_seq,
                      ticktext = paste0(100*beta_seq, "%"),
                      gridcolor = 'lightgray'),
         yaxis = list(title = "VaR",
                      gridcolor = 'lightgray'),
         plot_bgcolor = 'white',
         paper_bgcolor = 'white',
         hovermode = "x unified")

# Wykres 2: Oczekiwany Niedobór (95%)
es_plot <- plot_ly(results, x = ~beta, y = ~ES,
                   type = 'scatter', mode = 'lines+markers',
                   line = list(color = 'darkred', width = 3),
                   marker = list(size = 8, color = 'darkred'),
                   name = 'ES',
                   hovertemplate = paste(
                     'Beta: %{x:.1f}<br>',
                     'ES: %{y:.4f}<extra></extra>'
                   )) %>%
  layout(title = list(text = "Oczekiwany Niedobór (ES 95%)", x = 0.5),
         xaxis = list(title = "Waga na NVDA (β)",
                      tickmode = 'array',
                      tickvals = beta_seq,
                      ticktext = paste0(100*beta_seq, "%"),
                      gridcolor = 'lightgray'),
         yaxis = list(title = "ES",
                      gridcolor = 'lightgray'),
         plot_bgcolor = 'white',
         paper_bgcolor = 'white',
         hovermode = "x unified")

# OPCJA 1: Wyświetlanie wykresów obok siebie używając subplot
combined_plot <- subplot(var_plot, es_plot, 
                         nrows = 1, 
                         titleX = TRUE, 
                         titleY = TRUE,
                         margin = 0.05) %>%
  layout(title = list(text = "Miary Ryzyka Portfela vs Waga NVDA", 
                      x = 0.5, y = 0.98),
         showlegend = FALSE)

# Wyświetlenie połączonego wykresu
combined_plot
# Wyświetlenie tabeli wyników w lepszym formacie
# Tworzenie interaktywnej tabeli
datatable(results,
          colnames = c('Waga na NVDA (β)', 'VaR (95%)', 'ES (95%)'),
          options = list(
            pageLength = 11,  # Pokaż wszystkie wiersze
            dom = 't',        # Tylko pokaż tabelę
            columnDefs = list(
              list(className = 'dt-center', targets = 0:2),
              list(targets = 1:2, 
                   render = JS(
                     "function(data, type, row, meta) {
                       return type === 'display' ? 
                       parseFloat(data).toFixed(4) : data;
                     }"
                   ))
            )
          )) %>%
  formatPercentage('beta', 0) %>%
  formatRound(c('VaR', 'ES'), 4)