# SIMULACION AHORRO, EDAD DE PENSION Y TASAS DE INTERES
library(fitdistrplus)
dinero_a_gastar = 1.8*runif(1)
tiempo_gasto = 2

#tasa a la que empiezas a gastar el ahorro que se tiene
tasa_de_gasto = dinero_a_gastar/tiempo_gasto
# Paso de tiempo correspondiente a 1 mes en escala de años
mu = 1/12

# Se inicializa para conocer que personas alcanzan pensión
recibe_pension = 0
no_recibe_pension = 0
# Dinero disponible en el tiempo t. Por defecto esta vacio
Nt = c()

# Se analiza el número de pensionados, con base en años extra de vida
analisis_edad_pension <- function(extra_years){
  if (extra_years > 0){
    recibe_pension <<- recibe_pension + 1
  }else{
    no_recibe_pension <<- no_recibe_pension + 1
  }
}

# Ecuación en diferencias
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
  # Valor de vida que lanza la distribución Normal
  tiempo_vida = rnorm(1, mean = 63.977895, sd = 4.470362 )
  # tiempo extra de vida después de pensionarse
  extra_years = (tiempo_vida - edad_pension)
  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)*mu
  ahorro = 100
  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

    # Mientras f sea menos que u y tiempo sea menor que tiempo_vida
    while (F<U & tiempo<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))


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


# Datos en distribución logística
f2=fitdist(Nt,"logis")
summary(f2)
denscomp(f2, main = 'Logística')
qqcomp(list(f2), main = 'Logística')

gofstat(list(f1,f2)) #bondad del ajuste
qqcomp(list(f1,f2))

print(recibe_pension)
print(no_recibe_pension)