Tarea 3; Simulación

library(data.table)
library(lattice)
library(ggplot2)
library(caret)
library(jtools)
library(scales)
library(lattice)
library(MASS)
library(corrplot)

## corrplot 0.92 loaded

library(SimCorMultRes)

Pregunta I

Sección a)

# Función
f <- function(x) {
  ifelse(x > 0 & x <= 1, (3/5) * (4*x - x^2), 0)
}

# Data
data_df <- data.frame(x = seq(0, 1, by = 0.01),
                      y = f(seq(0, 1, by = 0.01)))

# Gráfico
ggplot() +
geom_line(data = data_df, aes(x = x, y = y), color = "blue") +
xlim(-1, 2) +
geom_segment(aes(x=-1,xend=0,y=0,yend=0),color='blue')+
 geom_segment(aes(x=1,xend=2,y=0,yend=0),color='blue')

Respuesta:

En esta figura se presenta la grafica asociada a la función f(x). En este grafico se considera la muestra de la función positiva y la sección igualada a cero. La forma en la que se presenta la grafica anterior se justifica en la función de densidad que se entrega por enunciado y fue comparada con la versión generada en otro software para justificar su intersección.

Adicionalmente, en este proceso se desarrolló una modificación al planteo general, para lograr que la función se mantuviera dentro de los lineamientos deseados, pero se considera que la función no considera el 1 dentro de su área distinta de 0.

Sección b)

# Planteo función t(x) as 3*x
t <- function(x) {
  return(3 * x)
}

# Planteo de integral c
c <- integrate(t, lower = 0, upper = 1)$value

# Planteo de r(x) as t(x)/c
r <- function(x) {
  return((3 * x)/c)
}

# Planteo de función f desde pregunta anterior
f <- function(x) {
  ifelse(x > 0 & x <= 1, (3/5) * (4*x - x^2), 0)
}

# Resulución de verificación para f(x)<=t(x)
x_values <- seq(0, 1, by = 0.01)
differences <- f(x_values)-t(x_values)

resultado=all(differences <= 0)


if (resultado == TRUE) {
  cat("Por lo tanto, según el resultado anterior se puede concluir que para todo x que pertenece a (0,1), f(x) <= t(x).\n")
} else {
  cat("Por lo tanto, según el resultado anterior, no se puede concluir que para todo x que pertenece a (0,1), f(x) <= t(x).\n")
}

## Por lo tanto, según el resultado anterior se puede concluir que para todo x que pertenece a (0,1), f(x) <= t(x).

# Funciones
f <- function(x) {
  ifelse(x > 0 & x <= 1, (3/5) * (4*x - x^2), 0)
}
t <- function(x) {
  return(3 * x)
}

# Data
data_df <- data.frame(x = seq(0, 1, by = 0.01),
                      y = f(seq(0, 1, by = 0.01)))
data_df1 <- data.frame(x = seq(0, 1, by = 0.01),
                      y = t(seq(0, 1, by = 0.01)))

# Gráfico
ggplot() +
geom_line(data = data_df, aes(x = x, y = y), color = "blue") +
  geom_line(data = data_df1, aes(x = x, y = y), color = "red")+
xlim(-1, 2) +
geom_segment(aes(x=-1,xend=0,y=0,yend=0),color='blue')+
 geom_segment(aes(x=1,xend=2,y=0,yend=0),color='blue') +   
geom_text(data = data.frame(x = c(1, 1), y = c(3, 1.8), label = c("t(x)", "f(x)")),
            aes(x = x, y = y, label = label), color = "black", size = 4)

Respuesta:

Dado el código anterior y la generación de esta nueva función, se puede determinar que la función f(x) es menor que la función nueva 3*x. Esto se justifica desde el código buscando el signo de la diferencia entre las funciones, pero también se justifica de manera grafica con la denotación de los gráficos en conjunto. En este caso, solo nos centraremos en la sección de x perteneciente a [1,0].

Sección c)

# Función de densidad de probabilidad g(x) = 1 en [0, 1]
g <- function(x) {
  ifelse(x >= 0 & x <= 1, 1, 0)
}

# Número de muestras deseadas
n_muestras <- 1000

# Inicializar un vector para almacenar las muestras aceptadas
muestras_aceptadas <- numeric(0)

# Contador de muestras aceptadas
muestras_aceptadas_count <- 0

# Realizar el método de aceptación/rechazo
while (muestras_aceptadas_count < n_muestras) {
  # Generar una muestra aleatoria de x en [0, 1]
  x <- runif(1)
  
  # Generar una muestra aleatoria de u en [0, 1]
  u <- runif(1)
  
  # Calcular f(x)
  fx <- f(x)
  
  # Calcular g(x)
  gx <- g(x)
  
  # Aceptar o rechazar la muestra según u <= f(x) / (c * g(x))
  if (u <= fx / (c*gx)) {
    muestras_aceptadas_count <- muestras_aceptadas_count + 1
    muestras_aceptadas <- c(muestras_aceptadas, x)
  }
}

# Calcular el ratio entre número de muestras aceptadas y número total de muestras
ratio_aceptadas_total <- length(muestras_aceptadas) / n_muestras

# Graficar el histograma de las muestras aceptadas
hist(muestras_aceptadas, breaks = 30, main = "Histograma de Muestras Aceptadas", xlab = "x")

# Imprimir el resultado

cat("Número de muestras necesarias para alcanzar 1000 muestras aceptadas:", length(muestras_aceptadas), "\n
")

## Número de muestras necesarias para alcanzar 1000 muestras aceptadas: 1000

cat("Ratio entre número de muestras aceptadas y número total de muestras:", ratio_aceptadas_total, "\n
")

## Ratio entre número de muestras aceptadas y número total de muestras: 1

cat("Número de muestras necesarias para alcanzar 1000 muestras aceptadas:", length(muestras_aceptadas), "\n
")

## Número de muestras necesarias para alcanzar 1000 muestras aceptadas: 1000

cat("Valor de c:", c, "\n")

## Valor de c: 1.5

Respuesta:

Dadas las condiciones denotadas en el código, y dado que no se genera una eliminación de observaciones, la cantidad de muestras requeridas para alcanzar las 1000 muestras aceptadas corresponde a 1000. Esto se justifica dado que la ratio de las muestras aceptadas es igual a 1 y el numero de muestras aceptadas es igual al numero de muestras generadas. En este caso se muestra que el método de aceptación generado puede ser reemplazado para mejorar la fidelidad de esta aceptación. No obstante, en este caso y resolución se opta por mantener el modelo aplicado para poder responder la pregunta, dado que se utilizan ambas funciones generadas en la pregunta anterior.

Sección d)

# Calcular f(x) simulada
x_values <- seq(0, 1, length.out = 1000)  # Valores de x en el rango [0, 1]

Fx <- numeric(length(x_values))

for (i in 1:length(x_values)) {
  Fx[i] <- integrate(f, 0, x_values[i])$value
}

# Calcular Fn(x) (distribución empírica acumulada)
Fn <- numeric(length(x_values))

for (i in 1:length(x_values)) {
  Fn[i] <- sum(muestras_aceptadas <= x_values[i]) / n_muestras
}

# Graficar F(x) y Fn(x) en el mismo gráfico
plot(x_values, Fx, type = "l", col = "blue", xlab = "x", ylab = "F(x)", main = "Distribución Acumulada")
lines(x_values, Fn, type = "l", col = "red")
legend("topright", legend = c("Dist. Simulada", "Dist. Empirica"), col = c("blue", "red"), lty = 1)

Respuesta:

Se puede observar que, dado el nivel de aceptación justificado anteriormente, las funciones de distribución acumulada entre la versión empírica y la simulada comparten varias características, dada la construcción de la simulación. En este caso, al no existir un nivel de rechazo tal de la muestra, se genera una distribución simulada similar a la versión empírica. Esta simulación representa la distribución empírica de una forma cercana, ya que buscamos generar una distribución que se acerque a las características planteadas.

Pregunta II

Sección a)

# Parámetros del modelo GBM
mu <- 0.15  # Tasa de crecimiento medio anual
sigma <- 0.30  # Volatilidad anual
P0 <- 50  # Precio inicial de una acción
n_years <- 30  # Número de años
n_simulations <- 1000  # Número de simulaciones
r <- 0.03  # Tasa de interés anual

# Función objetivo para simulación
simulate_portfolio <- function() {
  prices <- numeric(n_years + 1)
  prices[1] <- P0
  for (t in 1:n_years) {
    T <- t - 1  
    Z <- rnorm(1)
    investment <- 1000 * exp(r * T)
    prices[t + 1] <- (prices[t] + investment) * exp((mu - (sigma^2) / 2) + sigma * Z)
  }
  return(prices[n_years + 1])
}

# Implementación de simulación
set.seed(1234)
portfolio_values <- replicate(n_simulations, simulate_portfolio())
percentiles <- quantile(portfolio_values, c(0.05, 0.95))

# Recopilación de resultados
cat("Intervalo de predicción del 90% para el valor del fondo después de 30 años:\n
")

## Intervalo de predicción del 90% para el valor del fondo después de 30 años:

cat("Valor mínimo:", percentiles[1], "\n")

## Valor mínimo: 56958.9

cat("Valor máximo:", percentiles[2], "\n")

## Valor máximo: 2576016

Respuesta:

Los rangos presentados en este intervalo de predicción muestran el valor mínimo de la inversión junto al valor máximo luego de 30 años bajo las condiciones del enunciado. Si bien el riesgo asociado a la inversión es alto, se puede apreciar un aumento significativo en los montos, dada la constancia de la inversión. Este fondo, dado el riesgo asociado presenta una gran intervalo asociado a su rendimiento general.

Sección b)

mean_portfolio_value <- mean(portfolio_values)

se <- sd(portfolio_values) / sqrt(n_simulations)

t_critical <- qt(0.975, df = n_simulations - 1)

margin_of_error <- t_critical * se

lower_limit <- mean_portfolio_value - margin_of_error
upper_limit <- mean_portfolio_value + margin_of_error

half_interval <- (upper_limit - lower_limit) / 2
max_percent_deviation <- (half_interval / mean_portfolio_value) * 100

cat("Intervalo de confianza del 95% para el valor del fondo después de 30 años:\n
")

## Intervalo de confianza del 95% para el valor del fondo después de 30 años:

cat("Valor mínimo:",lower_limit , "\n")

## Valor mínimo: 648757

cat("Valor máximo:", upper_limit, "\n")

## Valor máximo: 810234.6

Respuesta:

Este intervalo de confianza acota más aun en rango de los valores, lo cual limita el máximo y mínimo a las condiciones determinadas de no más de un 0,5% del valor estimado. Por lo cual, es comprensible que se denote una disminución en el rango presentado en esta ocasión.

Sección c)

# Valor futuro del monto si no se invierte
PV_no_inversion <- 30000 * (1 + r)^30

# Simulación de los valores del fondo después de 30 años
simulated_values <- replicate(n_simulations, simulate_portfolio())

# Contar cuántas veces el valor del fondo es menor que el valor sin inversión
loss_count <- sum(simulated_values < PV_no_inversion)

# Calcular la probabilidad de pérdidas
probability_of_loss <- loss_count / n_simulations

# Imprimir la probabilidad de pérdidas
cat("Probabilidad de pérdidas después de 30 años:", probability_of_loss*100,"%")

## Probabilidad de pérdidas después de 30 años: 9.1 %

# ___________________________________________________

# Number of bootstrap samples
n_bootstrap_samples <- 1000

# Initialize a vector to store bootstrapped probabilities
bootstrapped_probabilities <- numeric(n_bootstrap_samples)

# Perform bootstrapping to estimate confidence interval
for (i in 1:n_bootstrap_samples) {
  # Resample with replacement from the simulated values
  resampled_values <- sample(simulated_values, replace = TRUE)
  
  # Count how many resampled values are less than PV_no_inversion
  bootstrapped_loss_count <- sum(resampled_values < PV_no_inversion)
  
  # Calculate the probability of loss for this resampled data
  bootstrapped_probability <- bootstrapped_loss_count / n_simulations
  
  # Store the bootstrapped probability
  bootstrapped_probabilities[i] <- bootstrapped_probability
}

# Calculate the lower and upper percentiles of the bootstrapped probabilities for the confidence interval
lower_percentile <- quantile(bootstrapped_probabilities, 0.025)
upper_percentile <- quantile(bootstrapped_probabilities, 0.975)

# Print the confidence interval
cat("95% Confidence Interval for Probability of Loss: [", lower_percentile, ", ", upper_percentile, "]\n")

## 95% Confidence Interval for Probability of Loss: [ 0.073 ,  0.108 ]

Respuesta:

Dado el volumen de la inversión y dada las condiciones del fondo se puede determinar que la probabilidad de perdidad frente a la volatibilidad del fonda A, esta en el rango de [ 0.073 , 0.108 ]. Este rango es reducido, lo cual puede denotar que existe una portafolio altamente diversificado el cual permite ir reduciendo la perdida dada la constante inversión de fondos.

Sección d)

Versión 2.a)

# Parámetros del modelo GBM
mu <- 0.15  # Tasa de crecimiento medio anual
sigma <- 0.30  # Volatilidad anual
P0 <- 50  # Precio inicial de una acción
n_years <- 30  # Número de años
n_simulations <- 1000  # Número de simulaciones
r <- 0.03  # Tasa de interés anual

# Función objetivo para simulación
simulate_portfolio <- function() {
  prices <- numeric(n_years + 1)
  prices[1] <- P0
  for (t in 1:n_years) {
    T <- t - 1  
    Z <- rnorm(1)
    investment <- 0 * exp(r * T)
    prices[t + 1] <- (prices[t] + investment) * exp((mu - (sigma^2) / 2) + sigma * Z)
  }
  return(prices[n_years + 1])
}

# Implementación de simulación
set.seed(1234)
portfolio_values <- replicate(n_simulations, simulate_portfolio())
percentiles <- quantile(portfolio_values, c(0.05, 0.95))

# Recopilación de resultados
cat("Intervalo de predicción del 90% para el valor del fondo después de 30 años:\n
")

## Intervalo de predicción del 90% para el valor del fondo después de 30 años:

cat("Valor mínimo:", percentiles[1], "\n")

## Valor mínimo: 88.24655

cat("Valor máximo:", percentiles[2], "\n")

## Valor máximo: 15744.7

Versión 2.b)

mean_portfolio_value <- mean(portfolio_values)

se <- sd(portfolio_values) / sqrt(n_simulations)

t_critical <- qt(0.975, df = n_simulations - 1)

margin_of_error <- t_critical * se

lower_limit <- mean_portfolio_value - margin_of_error
upper_limit <- mean_portfolio_value + margin_of_error

half_interval <- (upper_limit - lower_limit) / 2
max_percent_deviation <- (half_interval / mean_portfolio_value) * 100

cat("Intervalo de confianza del 95% para el valor del fondo después de 30 años:\n
")

## Intervalo de confianza del 95% para el valor del fondo después de 30 años:

cat("Valor mínimo:",lower_limit , "\n")

## Valor mínimo: 3388.257

cat("Valor máximo:", upper_limit, "\n")

## Valor máximo: 4462.473

Versión 2.c)

# Valor futuro del monto si no se invierte
PV_no_inversion <- 30000 * (1 + r)^30

# Simulación de los valores del fondo después de 30 años
simulated_values <- replicate(n_simulations, simulate_portfolio())

# Contar cuántas veces el valor del fondo es menor que el valor sin inversión
loss_count <- sum(simulated_values < PV_no_inversion)

# Calcular la probabilidad de pérdidas
probability_of_loss <- loss_count / n_simulations

# Imprimir la probabilidad de pérdidas
cat("Probabilidad de pérdidas después de 30 años:", probability_of_loss*100,"%")

## Probabilidad de pérdidas después de 30 años: 99.5 %

# ___________________________________________________

# Number of bootstrap samples
n_bootstrap_samples <- 1000

# Initialize a vector to store bootstrapped probabilities
bootstrapped_probabilities <- numeric(n_bootstrap_samples)

# Perform bootstrapping to estimate confidence interval
for (i in 1:n_bootstrap_samples) {
  # Resample with replacement from the simulated values
  resampled_values <- sample(simulated_values, replace = TRUE)
  
  # Count how many resampled values are less than PV_no_inversion
  bootstrapped_loss_count <- sum(resampled_values < PV_no_inversion)
  
  # Calculate the probability of loss for this resampled data
  bootstrapped_probability <- bootstrapped_loss_count / n_simulations
  
  # Store the bootstrapped probability
  bootstrapped_probabilities[i] <- bootstrapped_probability
}

# Calculate the lower and upper percentiles of the bootstrapped probabilities for the confidence interval
lower_percentile <- quantile(bootstrapped_probabilities, 0.025)
upper_percentile <- quantile(bootstrapped_probabilities, 0.975)

# Print the confidence interval
cat("95% Confidence Interval for Probability of Loss: [", lower_percentile, ", ", upper_percentile, "]\n")

## 95% Confidence Interval for Probability of Loss: [ 0.99 ,  0.999 ]

Respuesta:

En esta sección se presenta una reformulación de las preguntas anteriores bajo los cambios denotados por el fondo B. Esto permite generar una versión 2 de cada una de preguntas anteriores, ajustando los parámetros necesarios para el análisis correcto de cada uno de los fondos.

Sección e)

Respuesta:

Luego de realizar el ajuste pertinente al código para simular un Fondo B, podemos determinar que la inversión del dinero sin generar nuevo aporte al fondo va a aumentar la probabilidad de presentar perdida luego de 30 años. Esto permite que la elección frente a un Fondo sea clara dado el nivel de retorno presentado en el Fondo A y dada la probabilidad de perdida que se genera en el fondo B.

De esta forma, en mi opinión, el fondo A dado el análisis previo se ajusta de mejor manera a un fondo recomendable para la inversión.

Pregunta III

Reformulación Tarea 2

# Iniciando ejercicio :)
phi<-0.6
set.seed(1234)
N0 <- sample(5000, 1, replace = TRUE)


# Initialize variables for the simulation
N <- N0  # Assuming N is the same as N0 in your code
p_gorro_vector <- numeric(N)

# Matriz de correlación.
cor_matrix <- toeplitz(c(1, rep(phi, 2)))


for (i in 1:N0) {
  X1<- rexp(N0, 1)
  X2 <- rexp(N0, 1/2)
  X3 <- rnorm(N0, mean = 3, sd = 1)

    df<-data.table(X1,X2,X3)
    
  # Calculo de W
  W <- X1 + X2 + X3
  
  # Crear un vector Z para almacenar los resultados
  Z <- rep(0, N0)
  
  for (j in 1:N0) {
    if (W[j] > 10)
      Z[j] <- 1
    else
      Z[j] <- 0
  }
  
  # Calculate p_hat_g for this iteration
  p_gorro <- sum(Z) / N0
  
  # Store p_hat_g in the vector
  p_gorro_vector[i] <- p_gorro
}

# Calculate the mean and standard error
mean_p_gorro <- mean(p_gorro_vector)
se_p_gorro <- sqrt(mean(p_gorro_vector * (1 - p_gorro_vector)) / N)

# Calculate the 95% confidence interval
z_score <- qnorm(0.975)
lower_bound <- mean_p_gorro - z_score * se_p_gorro
upper_bound <- mean_p_gorro + z_score * se_p_gorro

# Print results
cat("Estimated probability p_hat_g:", mean_p_gorro, "\n")

## Estimated probability p_hat_g: 0.06701878

cat("Standard Error:", se_p_gorro, "\n")

## Standard Error: 0.007887995

cat("95% Confidence Interval: [", lower_bound, ",", upper_bound, "]\n")

## 95% Confidence Interval: [ 0.05155859 , 0.08247896 ]

Respuesta:

Esta reformulación del problema permite regenerar las variables necesarias para el análisis cuantitativo y presentar los resultados ajustados a la simulación original.

Sección a)

# Implementación de matriz de correlación
corrplot(cor_matrix, type="upper", order="hclust", sig.level = 0.01, insig = "blank")

Respuesta:

Esta matriz presenta la solicitud de la matriz de correlación, donde se asume que el nombre de C va a ser idéntico al de corr_matrix. Esta formulación presenta que las variables en la diagonal poseen una correlación de 100% dado que se refieren a las mismas variables y de 0.6 en el resto de la matriz.

Sección b)

set.seed(1234)
sample_size <- 5000
latent_correlation_matrix <- cor_matrix

# Formulación de distribuciones
different_marginal_distributions <- c("qexp", "qexp", "qnorm")

# Condiciónes de distribución original
qpars <- list(c(rate=1), 
              c(rate = 1/2),
              c(mean = 3, sd = 1))

# Aplicación de Rnorta
simulated <- rnorta(R = sample_size,
                    cor.matrix = latent_correlation_matrix,
                    distr = different_marginal_distributions,
                    qparameters = qpars)


# Resultados
cat("Matriz correlación simulada\n")

## Matriz correlación simulada

cor(simulated)

##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.5761787 0.5475926
## [2,] 0.5761787 1.0000000 0.5478259
## [3,] 0.5475926 0.5478259 1.0000000

cat("\nValores medios\n")

## 
## Valores medios

colMeans(simulated)

## [1] 0.9889279 2.0232557 3.0048725

cat("\nDesviación estándar\n")

## 
## Desviación estándar

apply(simulated, 2, sd)

## [1] 0.9755091 2.0243878 1.0029223

Respuesta:

Como podemos presentar, en la matriz de correlación simulada, los valores simulados son bastante cercanos a la solicitud. En este caso vemos una pequeña variación respecto a la correlación entre las variables y pequeñas variaciones entre los valores medios y la desviación estándar.

sección c)

# Calculate probability P(W > 10) and confidence interval for correlated samples
p_correlated <- mean(simulated[, 1] + simulated[, 2] + simulated[, 3] > 10)
conf_interval_correlated <- binom.test(sum(simulated[, 1] + simulated[, 2] + simulated[, 3] > 10), sample_size)$conf.int

# Calculate probability P(W > 10) and confidence interval for independent samples
p_independent <- mean(X1 + X2 + X3 > 10)
conf_interval_independent <- binom.test(sum(X1 + X2 + X3 > 10), sample_size)$conf.int

# Compare results
cat("Results for Correlated Samples:\n")

## Results for Correlated Samples:

cat("Probability P(W > 10):", p_correlated, "\n")

## Probability P(W > 10): 0.1192

cat("95% Confidence Interval:", conf_interval_correlated, "\n")

## 95% Confidence Interval: 0.1103429 0.1285043

cat("\nResults for Independent Samples:\n")

## 
## Results for Independent Samples:

cat("Probability P(W > 10):", p_independent, "\n")

## Probability P(W > 10): 0.07669323

cat("95% Confidence Interval:", conf_interval_independent, "\n")

## 95% Confidence Interval: 0.01217205 0.01921024

# Comment on the comparison

cat("Correlated vs. Independent Probability Difference:", p_correlated - p_independent, "\n")

## Correlated vs. Independent Probability Difference: 0.04250677

cat("Correlated vs. Independent Confidence Interval Difference:", conf_interval_correlated - conf_interval_independent, "\n")

## Correlated vs. Independent Confidence Interval Difference: 0.09817087 0.1092941

Respuesta:

Considerando el caso inicial de la formulación original, se puede realizar una comparación entre ambas instancias. En este caso, se utiliza una sample size de 5000 para la comparación. La comparación presenta una diferencia en la probabilidad, ya que en la muestra correlacionada es más probable que encontremos un valor de W superior a 10. Por otro lado, la longitud del intervalo de confianza es marginalmente superior en la muestra independiente, lo cual denota que hoy una ligera certeza mayor en la muestra correlacionada, pero esta afirmación debe ser justificada de forma estadística. Por lo tanto, existen ciertas diferencias al momento de generar las muestras y esto genera diferencias en los resultados, donde podemos presentar que el comportamiento correlacionado de la muestra entrega valores más cercanos a la realidad al mismo tiempo que se comporta con los requerimientos solicitados del enunciado.