# Required Libraries
library(copula)
library(MASS)
library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:MASS':
## 
##     select
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout

First: Define the marginal survival distributions (Weibull Distributed survival times)

Let X1 denote the RV for life 1 such that it is Weibull(3, 64) distributed Let X2 denote the RV for life 2 such that it is Weibull(8.2, 73) ditributed

# Set seed for reproducibility
set.seed(2005090192)

shape1 <- 3  # Shape parameter for Life 1
scale1 <- 64  # Scale parameter for Life 1

shape2 <- 8.2  # Shape parameter for Life 2
scale2 <- 73   # Scale parameter for Life 2

# Next: Simulate survival times for the marginal distributions
n <- 10000  # Number of samples
x1 <- rweibull(n, shape1, scale = scale1)
x2 <- rweibull(n, shape2, scale = scale2)

# Combine them into a data frame
data <- data.frame(Life1_Survival_Time = x1, Life2_Survival_Time = x2)

# A scatterplot of joint life survival times (plots X1 against X2)
ggplot(data, aes(x = Life1_Survival_Time, y = Life2_Survival_Time)) +
  geom_point(alpha = 0.7) +
  ggtitle("Scatterplot of Life Survival Times") +
  theme_minimal()

Second: Fit a Gumbel Copula to the data

gumbel_cop <- gumbelCopula(2)  # 2 is the default dependency parameter for the function

# Fit marginal distributions (Weibull) and estimate copula parameters
marginals <- mvdc(copula = gumbel_cop,
                  margins = c("weibull", "weibull"),
                  paramMargins = list(list(shape = shape1, scale = scale1),
                                      list(shape = shape2, scale = scale2)))

# Generate random sample pairs (u,v) from the copula 
## Note: this function obtains the pairs AND uses inversion to find the unstandardized values
samples <- rMvdc(n, marginals)

# Add the copula samples to a data frame
sample_data <- data.frame(Life1_Survival_Time = samples[, 1],
                          Life2_Survival_Time = samples[, 2])

# Visualize the accuracy of the estimated joint density (Using contour plots) from the Copula-Generated survival distribution times
ggplot(sample_data, aes(x = Life1_Survival_Time, y = Life2_Survival_Time)) +
  geom_point(alpha = 0.7) +
  geom_density2d() +
  ggtitle("Scatterplot of Life Survival Times with KDE Overlaid") +
  theme_minimal()

Third: Assess model fits for univariate distributions

# QQ plot for Life 1
## It is linear, so likely a good fit
qqplot(x1, rweibull(n, shape1, scale1), main = "QQ Plot for Life 1")
abline(0, 1, col = "red")

# QQ plot for Life 2
## It is linear, so likely a good fit
qqplot(x2, rweibull(n, shape2, scale2), main = "QQ Plot for Life 2")
abline(0, 1, col = "red")

Fourth: Visualize the joint distribution (as a 3D surface)

x_grid <- seq(min(x1), max(x1), length.out = 100)
y_grid <- seq(min(x2), max(x2), length.out = 100)

# Create meshgrid using expand.grid
cdf_grid <- expand.grid(x = x_grid, y = y_grid)

cdf_matrix <- as.matrix(cdf_grid)

# Calculate the CDF values for the Gumbel copula
cdf_values <- apply(cdf_matrix, 1, function(row) pMvdc(row, marginals))

# Reshape the CDF values into a grid for 3D visualization
Z <- matrix(cdf_values, nrow = length(x_grid), ncol = length(y_grid))

# Create a 3D surface plot of the CDF
plot_ly(x = ~x_grid, y = ~y_grid, z = ~Z, type = "surface") %>%
  layout(scene = list(xaxis = list(title = 'Life 1 Survival Time'),
                      yaxis = list(title = 'Life 2 Survival Time'),
                      zaxis = list(title = 'Cumulative Probability'),
                      title = 'Estimated Cumulative Distribution Function (CDF)'))

Fifth: Create a joint life probability/annuity table

# Fifth: Create a joint life probability/annuity table
age_x <- 50
age_y <- 55
age_limit <- 100

# Initialize a probability table
prob_table <- data.frame()

while(age_x <= age_limit & age_y <= age_limit) {
  prob_x_survives <- mean(sample_data$Life1_Survival_Time > age_x)
  prob_y_survives <- mean(sample_data$Life2_Survival_Time > age_y)
  joint_survival_prob <- mean(sample_data$Life1_Survival_Time > age_x & sample_data$Life2_Survival_Time > age_y)
  at_least_one_survives_prob <- prob_x_survives + prob_y_survives - joint_survival_prob
  
  prob_table <- rbind(prob_table, data.frame(Age_x = age_x,
                                             P_Life1_survives = prob_x_survives,
                                             Age_y = age_y,
                                             P_Life2_survives = prob_y_survives,
                                             Both_Survive = joint_survival_prob,
                                             At_least_one_survives = at_least_one_survives_prob))
  
  age_x <- age_x + 1
  age_y <- age_y + 1
}

head(prob_table, 20)
# Plot the probability table
ggplot(prob_table, aes(x = Age_x)) +
  geom_line(aes(y = P_Life1_survives, color = 'Life1')) +
  geom_line(aes(y = P_Life2_survives, color = 'Life2')) +
  geom_line(aes(y=At_least_one_survives, color='One Survives')) +
  geom_line(aes(y = Both_Survive, color = 'Joint')) +  # Corrected column name
  ggtitle('Survival Probabilities') +
  ylab('Probability') +
  theme_minimal()

Step 6: Create a table for Present Value of jointly contingent Insurance/Annuity Factors

# Set interest rate
interest_rate <- 0.05
v <- 1 / (1 + interest_rate)
d <- 1 - v

# Probabilities*discount
prob_table$P_Life_1_greater_x_v <- prob_table$P_Life1_survives * v
prob_table$P_Life_2_greater_y_v <- prob_table$P_Life2_survives * v
prob_table$P_Joint_Life1_greater_x_AND_Life2_greater_y_v <- prob_table$Both_Survive * v

# Single Life Annuities (ä_x, A_x for Life 1 and Life 2)
prob_table$ä_x <- rev(cumsum(rev(prob_table$P_Life_1_greater_x_v)))  # Cumulative sum in reverse order
prob_table$A_x <- 1 - prob_table$ä_x * d

prob_table$ä_y <- rev(cumsum(rev(prob_table$P_Life_2_greater_y_v)))  # Cumulative sum in reverse order
prob_table$A_y <- 1 - prob_table$ä_y * d

# Joint Life Annuity (ä_(xy), A_(xy) for joint lives)
prob_table$ä_xy <- rev(cumsum(rev(prob_table$P_Joint_Life1_greater_x_AND_Life2_greater_y_v)))
prob_table$A_xy <- 1 - prob_table$ä_xy * d

# Last Survivor (ä_(xy:bar), A_(xy:bar))
prob_table$ä_xy_bar <- prob_table$ä_x + prob_table$ä_y - prob_table$ä_xy
prob_table$A_xy_bar <- prob_table$A_x + prob_table$A_y - prob_table$A_xy

# Reversionary Annuity (x survives y dies, ä_y|x)
prob_table$ä_y_given_x <- prob_table$ä_y - prob_table$ä_xy

# Reversionary Annuity (y survives x dies, ä_x|y)
prob_table$ä_x_given_y <- prob_table$ä_x - prob_table$ä_xy

# Display the first 20 rows of the result
head(prob_table, 20)
# Plotting the results (for each of the calculated columns)
library(ggplot2)

columns_to_plot <- c("ä_x", "A_x", "ä_y", "A_y", "ä_xy", "A_xy", "ä_xy_bar", "A_xy_bar", "ä_y_given_x", "ä_x_given_y")

for (column in columns_to_plot) {
  plot <- ggplot(prob_table, aes(x = 1:nrow(prob_table), y = get(column))) +
    geom_line() +
    ggtitle(paste("Plot of", column)) +
    xlab("Index") +
    ylab(column) +
    theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5))
  
  print(plot)  # Explicitly print the ggplot object
}