######################
## Simulazione in R ##
######################

# 1) caricare il file IPS55.TAV contenente valori di lx e simulare la durata di vita troncata (per eccesso o per difetto)
setwd("inserire qui la directory") # oppure in Rstudio -> Session ->
IPS55M <- read.table(file = "IPS55M.TAV", header = FALSE, col.names = c("x", "lx"))

tpx <- IPS55M$lx / IPS55M$lx[1] # prob di sopravvivenza
tpx <- c(tpx, 0) # aggiungi uno 0 alla fine
x <- c(IPS55M$x, tail(IPS55M$x, n = 1) + 1) # aggiungi l'eta' estrema (prima eta' che non e' possibile raggiungere in vita)



# RNGkind
# normal.kind permette di scegliere il metodo per simulare variabili aleatorie normali
# "Kinderman-Ramage",
# "Buggy Kinderman-Ramage"
# "Ahrens-Dieter"
# "Box-Muller"
# "Inversion" (default)
# "user-supplied"
# 2) verificare il tempo impiegato per 10 ^ 7 simulazioni da  Box-Muller e dal metodo di default



# ATTENZIONE ALLA PARAMETRIZZAZIONE DELLE DISTRIBUZIONI
# 3) Esempio: nella distribuzione gamma si deve specificare "shape" e uno tra "rate" o "scale"
# rate = 1 / scale


?rgamma

# 4) Scrivere una funzione per generare la somma di lognormali (iid)
# osservazione: non e' nota la distribuzione di tale variabile aleatoria





# 5) Calcolare E[X], E[X ^ 2], P(X > 8) tramite Monte Carlo
# quando X = somma di 5 lognormali(0, 1)
EX.true <- 5 * exp(0 + 0.5 * 1 ^ 2)


# stima, intervallo di confidenza

# ampiezza intervallo

# errore assoluto

# errore relativo

# svolgere il n 6) e poi continuare


# 6) scrivere funzione che dato un campione simulato, calcola il valor medio e l'intervallo di confidenza
MonteCarlo <- function(x, level = 0.95)
{
  nsim <- length(x)
  MC <- mean(x)
  SE <- sd(x) / sqrt(nsim)
  z <- qnorm((1 - level) / 2, lower.tail = FALSE)
  MC.l <- MC - z * SE
  MC.u <- MC + z * SE
  IC <- c(MC.l, MC.u)
  names(IC) <- as.character(c((1 - level) / 2, level + (1 - level) / 2))
  return(list(MC, IC))
}

nsim1 <- 10 ^ (1 : 6) # da 10 a 1 milione
XX <- matrix(nrow = 6, ncol = 6) # risultati
for (i in seq_along(nsim1))
{
  set.seed(1)
  X <- rSommaLognormali(nsim = nsim1[i], nl = 5, meanlog = 0, sdlog = 1)
  MC <- MonteCarlo(X, level = 0.99)
  XX[i, ] <- c(MC[[1]], MC[[2]], MC[[2]][2] - MC[[2]][1], abs(MC[[1]] - EX.true), 100 * abs((MC[[1]] - EX.true) / EX.true))
}

risultato <- data.frame(nsim = nsim1,
                        MC = XX[, 1],
                        inf = XX[, 2],
                        sup = XX[, 3],
                        delta = XX[, 4],
                        errore.a = XX[, 5],
                        errore.r = XX[, 6])

risultato

# 7) Scrivere funzione per la f. di ripartizione, densita', quantile e simulazione delle vita residua di un individuo distribuita come una Gompertz
# usare x = 40, alpha = 2.67433 * 10 ^ -5, beta = log(1.098)
# confrontare la distribuzione di Tx con density() vs densita' esatta
# calcolare valor medio di T con Monte Carlo e via integrazione numerica della funzione di sopravvivenza

pGompertz <- function(t, x, alpha, beta, lower.tail = TRUE)
{
  B <- exp(beta * x)
  p <- 1 - exp(- alpha * B * (exp(beta * t) - 1) / beta)
  if (lower.tail) return(p) else return(1 - p)
}

dGompertz <- function(t, x, alpha, beta)
{
  B <- exp(beta * x)
  return(alpha * B * exp(beta * t - alpha * B * (exp(beta * t) - 1) / beta))
}

qGompertz <- function(p, x, alpha, beta)
{
  B <- exp(beta * x)
  log(1 - beta * log(1 - p) / (alpha * B)) / beta
}






tt <- 0 : 80
plot(tt, dGompertz(t = tt, x = 40, alpha = 2.67433 * 10 ^ -5, beta = log(1.098)), type = "l", xlab = "t", ylab = "Densita'", main = "densita' della Gompertz")

hist(T, freq = FALSE, breaks = 100)
lines(density(T), col = "red")


# calcolo via Monte Carlo

# via integrazione numerica della funzione di sopravvivenza


# errore assoluto

# errore relativo


# per diminuire l'ampiezza dell'intervallo di confidenza di fattore 10, aumentiamo di 100 il numero di smiulazioni



# 8) calcolare valore attuale atteso di

# caso morte vita intera
# temporanea caso morte
# capitale differito
# mista
# rendita vitalizia (continua)

# usando il modello in 7) e tasso d'interesse 2%

set.seed(1)
nsim <- 10000
T <- rGompertz(nsim = nsim,
                     x = 40,
                     alpha = 2.67433 * 10 ^ -5,
                     beta = log(1.098))

i <- 0.02
v <- 1 / (1 + i)

# caso morte vita intera



# temporanea caso morte, scadenza 20 anni


# capitale differito, scadenza 10 anni


# mista


# rendita vitalizia (continua)




# 9) la legge forte dei grandi numeri (e Monte Carlo) non funziona con distribuzioni con media non definita
# provare a calcolare il valore atteso di una v.a. con distribuzione di Cauchy e confrontarlo con un'esponenziale




# 10) metodo accetazione rifiuto per simulare una durata di vita con intensita' Makeham

rMakeham <- function(nsim, x, alpha, beta, gamma)
{
  Tx <- vector(mode = "numeric", length = nsim)
  count <- 0
  B <- exp(beta * x)
  C <- 1 + gamma / (alpha * B)
  repeat
  {
    Z1 <- rGompertz(nsim = 1, x, alpha, beta)
    U <- runif(1)
    f.Cg <- (gamma * exp(- beta * Z1) + alpha * B)/(gamma + alpha * B)* exp(- gamma * Z1)
    if (U <= f.Cg)
    {
      count <- count + 1
      Tx[count] <- Z1
      if (count == nsim) break
    }
  }
  return(Tx)
}

x <- 40
alpha <- 2.67433 * 10 ^ -5
beta <- log(1.098)
gamma <- 0.001

debug(rMakeham)
TT <- rMakeham(nsim = 1000, x, alpha, beta, gamma)
MonteCarlo(TT)

dMakeham <- function(t, x, alpha, beta, gamma)
{
  B <- exp(beta * x)
  return((gamma + alpha * B * exp(beta * t)) * exp(- gamma * t - ((alpha * B) / beta) * (exp(beta * t) - 1)))
}

integrate(function(t)(t * dMakeham(t, x = x, alpha = alpha, beta = beta, gamma = gamma)), lower = 0, upper = 80)

# 11) Calcolo di un integrale via montecarlo
# Esempio: integrale di exp(-x ^ 2) / (x ^ 2 + 1)
# su intervalli limitati e illimitati


# integrazione su intervalli (limitati) diversi
FUNab <- function(x, a, b, FUN)(b - a) * FUN(a + (b - a) * x)

# su [1,2]


# su [-0.1,0.1]

# integrale su intervallo illimitato
# cambi di variabile:
# y = (exp(x) + 1) ^ -1
# oppure
# y = (arctg(x) - pi / 2) / pi



# 12) calcolo della funzione gamma
# Gamma(a) = integrale su [0,+Inf) di x ^ (a - 1) * exp(- x)
# cambio di variabile y = exp(- x)


# cambio di variabile y = 1 / (x + 1)


# 13) Calcolo di integrali multidimensionali (su ipercubi) via MC
# Esempio: f(x) = exp( - ||x|| ^ 2) integrato su [0,1] x ... x [0,1]
FUN <- function(x)exp(- x %*% x)
FUN1 <- function(x)exp(- sum(x ^ 2))

# valore corretto
INT <- function(n)(sqrt(pi) * (pnorm(sqrt(2)) - 0.5)) ^ n

# con il pacchetto "cubature"
# install.packages("cubature")
library(cubature)
library(help = cubature)
?hcubature
NUM <- function(n, ...)hcubature(f = FUN, lowerLimit = rep(0, n), upperLimit = rep(1, n), ...)

# via MC


# confronto

# dimensione 1

# dimensione 2

# dimensione 3

# dimensione 7

# dimensione 8

# dimensione 10

# dimensione 20