# SIMULACION AHORRO, EDAD DE PENSION Y TASAS DE INTERES
library(fitdistrplus)

dinero_a_gastar = runif(1, min = 1, max = 3)
tiempo_gasto = 2

tasa_de_gasto = dinero_a_gastar/tiempo_gasto #tasa a la que empiezas a gastar el ahorro que se tiene, por ejemplo 1 millon cada tres meses. ¿cuanto te demoras gastando un millon?

mu = 1/12 # Paso de tiempo correspondiente a 1 mes en escala de años

recibe_pension = 0
no_recibe_pension = 0

Nt = c() # Dinero disponible en el tiempo t. Por defecto esta vacio

analisis_edad_pension <- function(edad){
  if (edad>0){
    recibe_pension <<- recibe_pension + 1
  }else{
    no_recibe_pension <<- no_recibe_pension + 1
  }
}

differential_equation <- function(F, mu, tasa_de_gasto, ahorro){
  return(F+mu*(1-F)*tasa_de_gasto*ahorro)
}

for (iteracion in 1:1000)
{
  # TIEMPOS

  edad_pension = 62
  beneficiario = rlogis(1, location = 64.13, scale = 2.67)  # Tiempo que un beneficiario recibe pension de uno ya fallecido
  tiempo_vida = rlogis(1, location = 64.13, scale = 2.67) # Valor de vida que lanza la distribución logistica
  extra_years = (tiempo_vida + beneficiario - edad_pension)*mu # tiempo extra de vida después de pensionarse y en caso de muerte recibe pension un beneficiario
  tiempo = 0 # tiempo que va a transcurrir durante la iteración

  # VALORES QUE AFECTAN EL DINERO
  #interes = rnorm(1, mean = 0.02780, sd = 0.0127)
  #interes =
  interes = rlogis(1, location = 0.02628266, scale = 0.00669499)*mu

  ahorro = 100 # Arranca en 100 en todos los casos, se puede modificar

  analisis_edad_pension(extra_years)

  while (tiempo < extra_years & ahorro > 0) {

    U = runif(1) # Aleatorio entre 0 y 1 que se usa para comparar con F
    F = mu # Arranca en mu por que es el primer periodo de tiempo que se afecta

    while (F<U & tiempo<tiempo_vida) { #mientras f sea menos que u y tiempo sea menor que tiempo_vida
      F = differential_equation(F, mu, tasa_de_gasto, ahorro)
      tiempo = tiempo+mu
    }

    ahorro = (ahorro-dinero_a_gastar)*(1+interes)


  }
  Nt[iteracion] = ahorro
  print("Años extra de vida")
  print(extra_years)

}
Nt
hist(Nt, main = "Ahorro en Banco")
quantile(Nt,c(0.025,0.50,0.975))

normal = fitdist(Nt,"norm")
summary(normal)
denscomp(normal, main = "Compara ahorro con la distribución Normal")
qqcomp(list(normal), main = 'normal')

exponencial = fitdist(Nt,"exp")
summary(exponencial)
denscomp(exponencial, main = "Compara ahorro con la distribución Exponencial")
qqcomp(list(exponencial), main = 'Exponencial')

lnormal = fitdist(Nt,"lnorm")
summary(lnormal)
denscomp(lnormal, main = "Compara ahorro con la distribución lnormal")
qqcomp(list(lnormal), main = 'lnormal')

#gofstat(list(normal, exponencial, lnormal)) #bondad del ajust


print(recibe_pension)
print(no_recibe_pension)