Math 563 Final Project Computations

1: PELVE Distribution Applications

1.1: Standard Normal Distribution Application

1.1.1: PELVE Function

normal_pelve = function(epsilon, mu, sigma_sq) {
  
  modified_d_norm = function(c_val) {
    temp_qnorm = qnorm(1 - c_val * epsilon)
    return(dnorm(temp_qnorm) / (c_val * epsilon))
  }
  
  # Assume epsilon <= 1
  c_grid = seq(from = 1, 
               to = 1 / epsilon, 
               length.out = 1000 * ceiling(1 / epsilon))
  c_solution = Inf
  
  VaR = qnorm(1 - epsilon)
  
  ES_c_grid = lapply(c_grid,
                     modified_d_norm)
  
  for (i in 1:length(ES_c_grid)) {
    if (ES_c_grid[[i]] <= VaR) {
      c_solution = min(c_solution, c_grid[i])
    }
  }
  
  return(c_solution)
}

1.1.2: Mean Testing

# Testing for standard normal
normal_pelve(0.1, 0, 1)

## [1] 2.457246

# Testing for translated normal (non-zero mean and variance 1)
means = seq(from = -1.5,
            to = 1,
            length.out = 200)
pelve_means = unlist(lapply(means,
                            epsilon = 0.1,
                            sigma_sq = 1,
                            normal_pelve))
plot(pelve_means ~ means,
     main = "A Plot of PELVE vs Mean for the Normal Distribution",
     xlab = "PELVE Value",
     ylab = "Mean Value")

1.1.3: Variance Testing

# Testing for scales normal (zero mean, non-1 variance
variance = seq(from = 0,
               to = 1.5,
               length.out = 200)
pelve_variances = unlist(lapply(variance,
                                epsilon = 0.1,
                                mu = 0,
                                normal_pelve))
plot(pelve_variances ~ variance,
     main = "A Plot of PELVE vs Variance for the Normal Distribution",
     xlab = "PELVE Value",
     ylab = "Variance Value")

PELVE peaks at a point; investigate what is going on here.

1.3: Student’s T-distribution Application

1.3.1: PELVE Function

t_pelve = function(epsilon, freedom) {
  
  ES = function(c_val) {
    first_half_ES = 1 / (c_val * epsilon) * dt(qt(1 - c_val * epsilon, 
                                                  df = freedom),
                                               df = freedom)
    second_half_ES = (freedom + qt(1 - c_val * epsilon,
                                   df = freedom) ^ 2) / (freedom - 1)
    
    return(first_half_ES * second_half_ES)
  }
  
  # Assume epsilon <= 1
  c_grid = seq(from = 1, 
               to = 1 / epsilon, 
               length.out = 1000 * ceiling(1 / epsilon))
  c_solution = Inf
  
  VaR = qt(1 - epsilon,
           df = freedom)
  
  ES_c_grid = lapply(c_grid,
                     ES)

  for (i in 1:length(ES_c_grid)) {
    if (is.na(ES_c_grid[[i]])) {
      break
    }
    if (ES_c_grid[[i]] <= VaR) {
      c_solution = min(c_solution, c_grid[i])
    }
  }
  
  return(c_solution)
}

1.3.2: Student’s t-distribution PELVE Paper Computation Validation

# freedom 2
epsilons = c(0.1, 0.05, 0.01, 0.005)
t_epsilons_2 = unlist(lapply(epsilons,
                             freedom = 2,
                             t_pelve))

t_data_frame_2 = data.frame(epsilons, t_epsilons_2)
print(t_data_frame_2)

##   epsilons t_epsilons_2
## 1    0.100      3.60036
## 2    0.050      3.80074
## 3    0.010      3.96013
## 4    0.005      3.98004

# freedom 10
t_epsilons_10 = unlist(lapply(epsilons,
                              freedom = 10,
                              t_pelve))

t_data_frame_10 = data.frame(epsilons, t_epsilons_10)
print(t_data_frame_10)

##   epsilons t_epsilons_10
## 1    0.100      2.582358
## 2    0.050      2.654983
## 3    0.010      2.745387
## 4    0.005      2.768124

# freedom 30
t_epsilons_30 = unlist(lapply(epsilons,
                              freedom = 30,
                              t_pelve))

t_data_frame_30 = data.frame(epsilons, t_epsilons_30)
print(t_data_frame_30)

##   epsilons t_epsilons_30
## 1    0.100      2.495050
## 2    0.050      2.554278
## 3    0.010      2.629556
## 4    0.005      2.648723

1.3.3: Error Testing

epsilons = seq(from = 0.0001,
               to = 0.5,
               length.out = 100)
pelve_error = unlist(lapply(epsilons,
                            freedom = 10,
                            t_pelve))
plot(pelve_error ~ epsilons,
     main = "A Plot of PELVE vs Error for the Student's t-distribution",
     xlab = "Epsilon Value",
     ylab = "PELVE Value")

1.4: Rayleigh Distribution Application

1.4.1: PELVE Function

rayleigh_pelve = function(epsilon) {
  
  erf = function(x) 2 * pnorm(x * sqrt(2)) - 1
  
  ES = function(c_val) {
    # Calculated explicit expression through integral-calculator.com
    first_numerator = sqrt(pi) * erf(sqrt(-log(c_val * epsilon)))
    second_numerator = 2 * (c_val * epsilon) * sqrt(-log(c_val * epsilon))
    return((first_numerator + second_numerator) / (sqrt(2) * c_val * epsilon))
  }
  
  # Assume epsilon <= 1
  c_grid = seq(from = 1, 
               to = 1 / epsilon, 
               length.out = 1000 * ceiling(1 / epsilon))
  c_solution = Inf
  
  VaR = sqrt(-2 * log(epsilon))
  
  ES_c_grid = lapply(c_grid,
                     ES)

  for (i in 1:length(ES_c_grid)) {
    if (is.na(ES_c_grid[[i]])) {
      break
    }
    if (ES_c_grid[[i]] <= VaR) {
      c_solution = min(c_solution, c_grid[i])
    }
  }
  
  return(c_solution)
}

1.4.2: Rayleigh Testing

# variance 0.5
epsilons = c(0.1, 0.05, 0.01, 0.005)
rayleigh_pelve = unlist(lapply(epsilons,
                                  rayleigh_pelve))

rayleigh_data_frame_small = data.frame(epsilons, rayleigh_pelve)
print(rayleigh_data_frame_small)

##   epsilons rayleigh_pelve
## 1    0.100       6.556256
## 2    0.050      11.999650
## 3    0.010      50.700477
## 4    0.005      95.431942