1º - Um modelo de múltiplo decremento com duas causas de decremento é especificado pelas seguintes forças de decrementos:

Decremento 1: \(μ{_x^1}(t) = \frac{1}{100-(x+t)}, t\le100-x\)
Decremento 2: \(μ{_x^2}(t) = \frac{2}{100-(x+t)}, t\le100-x\)

Indivíduo de idade x = 50 anos.

library(readxl)
library(DT)
library(knitr)

x <- 50
t <- seq(0, 49, 1)

a) Fdp da distribuição conjunta

\(f_{T, J}(t, j) = {_t}p{_x^{(\tau)}} \cdot \mu{_x^{(j)}}(t)\)

Forças decrementais:

Decremento 1:
\(\mu{_{x}^{1}} (t) = \frac{1}{100-(x+t)}\)
\(\mu{_{50}^{1}} (t) = \frac{1}{100-(50+t)} = \frac{1}{50-t}\)

μ1 <- 1 / (50 - t); μ1
##  [1] 0.02000000 0.02040816 0.02083333 0.02127660 0.02173913 0.02222222
##  [7] 0.02272727 0.02325581 0.02380952 0.02439024 0.02500000 0.02564103
## [13] 0.02631579 0.02702703 0.02777778 0.02857143 0.02941176 0.03030303
## [19] 0.03125000 0.03225806 0.03333333 0.03448276 0.03571429 0.03703704
## [25] 0.03846154 0.04000000 0.04166667 0.04347826 0.04545455 0.04761905
## [31] 0.05000000 0.05263158 0.05555556 0.05882353 0.06250000 0.06666667
## [37] 0.07142857 0.07692308 0.08333333 0.09090909 0.10000000 0.11111111
## [43] 0.12500000 0.14285714 0.16666667 0.20000000 0.25000000 0.33333333
## [49] 0.50000000 1.00000000

Decremento 2:
\(\mu{_{x}^{2}} (t) = \frac{2}{100-(x+t)}\)
\(\mu{_{x}^{2}} (t) = \frac{2}{100-(50+t)} = \frac{2}{50-t}\)

μ2 <- 2 / (50 - t); μ2
##  [1] 0.04000000 0.04081633 0.04166667 0.04255319 0.04347826 0.04444444
##  [7] 0.04545455 0.04651163 0.04761905 0.04878049 0.05000000 0.05128205
## [13] 0.05263158 0.05405405 0.05555556 0.05714286 0.05882353 0.06060606
## [19] 0.06250000 0.06451613 0.06666667 0.06896552 0.07142857 0.07407407
## [25] 0.07692308 0.08000000 0.08333333 0.08695652 0.09090909 0.09523810
## [31] 0.10000000 0.10526316 0.11111111 0.11764706 0.12500000 0.13333333
## [37] 0.14285714 0.15384615 0.16666667 0.18181818 0.20000000 0.22222222
## [43] 0.25000000 0.28571429 0.33333333 0.40000000 0.50000000 0.66666667
## [49] 1.00000000 2.00000000
Força total de decremento:

\(\mu{_{x}^{(\tau)}}(t) = \sum{_{j=1}^{m}}\mu{_{x}^{(j)}}(t)\)
\(\mu{_{50}^{(\tau)}}(t) = \mu{_{50}^{1}}(t) + \mu{_{50}^{2}}(t)\)
\(\mu{_{50}^{(\tau)}}(t) = \frac{1}{50-t} + \frac{2}{50-t}\)
\(\mu{_{50}^{(\tau)}}(t) = \frac{3}{50-t}\)

μτ <- μ1 + μ2; μτ
##  [1] 0.06000000 0.06122449 0.06250000 0.06382979 0.06521739 0.06666667
##  [7] 0.06818182 0.06976744 0.07142857 0.07317073 0.07500000 0.07692308
## [13] 0.07894737 0.08108108 0.08333333 0.08571429 0.08823529 0.09090909
## [19] 0.09375000 0.09677419 0.10000000 0.10344828 0.10714286 0.11111111
## [25] 0.11538462 0.12000000 0.12500000 0.13043478 0.13636364 0.14285714
## [31] 0.15000000 0.15789474 0.16666667 0.17647059 0.18750000 0.20000000
## [37] 0.21428571 0.23076923 0.25000000 0.27272727 0.30000000 0.33333333
## [43] 0.37500000 0.42857143 0.50000000 0.60000000 0.75000000 1.00000000
## [49] 1.50000000 3.00000000
μτ <- 3 / (50 - t); μτ
##  [1] 0.06000000 0.06122449 0.06250000 0.06382979 0.06521739 0.06666667
##  [7] 0.06818182 0.06976744 0.07142857 0.07317073 0.07500000 0.07692308
## [13] 0.07894737 0.08108108 0.08333333 0.08571429 0.08823529 0.09090909
## [19] 0.09375000 0.09677419 0.10000000 0.10344828 0.10714286 0.11111111
## [25] 0.11538462 0.12000000 0.12500000 0.13043478 0.13636364 0.14285714
## [31] 0.15000000 0.15789474 0.16666667 0.17647059 0.18750000 0.20000000
## [37] 0.21428571 0.23076923 0.25000000 0.27272727 0.30000000 0.33333333
## [43] 0.37500000 0.42857143 0.50000000 0.60000000 0.75000000 1.00000000
## [49] 1.50000000 3.00000000
Encontrando a probabilidade de que o indivíduo chegue até o tempo t sem sofrer nenhum decremento:

\({_t}p{_{x}^{(\tau)}} = e^{-\int_{0}^{t}\mu{_{x}^{(\tau)}}(t)dt}\)
\({_t}p{_{50}^{(\tau)}} = e^{-\int_{0}^{t}\mu{_{50}^{(\tau)}}(t)dt}\)
Substituindo a força total de decremento:
\({_t}p{_{50}^{(\tau)}} = e^{-\int_{0}^{t}\frac{3}{50 - t}dt}\)
Resolvendo a integral:
\(\int_{0}^{t}\frac{3}{50 - t} dt = - 3ln(50 - t) + 3ln(50)\)
Voltando a expressão da probabilidade:
\({_t}p{_{50}^{(\tau)}} = e^{-(-3ln(50 - t) + 3ln(50))}\)
\({_t}p{_{50}^{(\tau)}} = e^{3ln(50 - t) - 3ln(50)}\)
Aplicando a propriedade dos logaritmos, \(ln(a) - ln(b) = ln(\frac{a}{b})\):
\({_t}p{_{50}^{(\tau)}} = e^{3ln\frac{(50-t)}{50}}\)
\({_t}p{_{50}^{(\tau)}} = (\frac{(50-t)}{50})^3\)

pτ <- ((50 - t) / 50) ^ 3; pτ
##  [1] 1.000000 0.941192 0.884736 0.830584 0.778688 0.729000 0.681472 0.636056
##  [9] 0.592704 0.551368 0.512000 0.474552 0.438976 0.405224 0.373248 0.343000
## [17] 0.314432 0.287496 0.262144 0.238328 0.216000 0.195112 0.175616 0.157464
## [25] 0.140608 0.125000 0.110592 0.097336 0.085184 0.074088 0.064000 0.054872
## [33] 0.046656 0.039304 0.032768 0.027000 0.021952 0.017576 0.013824 0.010648
## [41] 0.008000 0.005832 0.004096 0.002744 0.001728 0.001000 0.000512 0.000216
## [49] 0.000064 0.000008

Expressões para a fdp da distribuição conjunta de T e J:

\(f_{T, J}(t, j) = {_t}p{_x^{(\tau)}} \cdot \mu{_x^{(j)}}(t)\)

\[ f_{T,J}(t,j) = \begin{cases} f_{T,J}(t,1) = (\frac{(50-t)}{50})^3 \cdot \frac{1}{50-t} = \frac{(50 - t)^2}{50^3}, t\le50 \\ f_{T,J}(t,2) = (\frac{(50-t)}{50})^3 \cdot \frac{2}{50-t} = \frac{2 \cdot (50 - t)^2}{50^3}, t\le50 \end{cases} \]

b) Fdp da distribuição marginal \(f_{T}(t)\)

\(f_{T}(t) = {_t}p{_x^{(\tau)}} \cdot \mu{_x^{(\tau)}}(t)\)
\(f_{T}(t) = (\frac{(50-t)}{50})^3 \cdot \frac{3}{50-t} = \frac{3(50-t)^2}{50^3}, t\le50\)

fT <- 3 * (50 - t) ^ 2 / 50 ^ 3; fT
##  [1] 0.060000 0.057624 0.055296 0.053016 0.050784 0.048600 0.046464 0.044376
##  [9] 0.042336 0.040344 0.038400 0.036504 0.034656 0.032856 0.031104 0.029400
## [17] 0.027744 0.026136 0.024576 0.023064 0.021600 0.020184 0.018816 0.017496
## [25] 0.016224 0.015000 0.013824 0.012696 0.011616 0.010584 0.009600 0.008664
## [33] 0.007776 0.006936 0.006144 0.005400 0.004704 0.004056 0.003456 0.002904
## [41] 0.002400 0.001944 0.001536 0.001176 0.000864 0.000600 0.000384 0.000216
## [49] 0.000096 0.000024
fT <- pτ * μτ; fT
##  [1] 0.060000 0.057624 0.055296 0.053016 0.050784 0.048600 0.046464 0.044376
##  [9] 0.042336 0.040344 0.038400 0.036504 0.034656 0.032856 0.031104 0.029400
## [17] 0.027744 0.026136 0.024576 0.023064 0.021600 0.020184 0.018816 0.017496
## [25] 0.016224 0.015000 0.013824 0.012696 0.011616 0.010584 0.009600 0.008664
## [33] 0.007776 0.006936 0.006144 0.005400 0.004704 0.004056 0.003456 0.002904
## [41] 0.002400 0.001944 0.001536 0.001176 0.000864 0.000600 0.000384 0.000216
## [49] 0.000096 0.000024

c) Fdp da distribuição marginal \(f_{J}(j)\)

\(f_{J}(j) = {_\infty}q{_{x}^{(j)}}\)
Calculando a probabilidade de que o decremento ocorra até o tempo t
\(_{50}q{_{50}^{(1)}} = \int{_{0}^{50}}\frac{(50 - t)^2}{50^3}dt = \frac{1}{3}\)

q1 <- function(t) {
  ((50 - t) ^ 2 / 50 ^ 3)
}
funcaoq1 <- integrate(q1, 0, 50)$value; funcaoq1
## [1] 0.3333333

\(_{50}q{_{50}^{(2)}} = 1 - {_{50}}q{_{50}^{(1)}} = \frac{2}{3}\)

q2 <- 1 - funcaoq1; q2
## [1] 0.6666667

d) Fdp da distribuição condicional de \(J\), dado um decremento no tempo \(t - f_{J|T}(j|t)\) 

\(f_{J|T}(j|t) = \frac{f_{T, J}(t, j)}{f_{T}(t)} = \frac{{_t}p{_x^{(\tau)}} \cdot \mu{_x^{(j)}}(t)}{{_t}p{_x^{(\tau)}} \cdot \mu{_x^{(\tau)}}(t)} = \frac{\mu_{x}^{(j)}(t)}{\mu_{x}^{(\tau)}(t)}\)
\(f_{J|T}(1|t) = \frac{\mu_{x}^{(1)}(t)}{\mu_{x}^{(\tau)}(t)} = \frac{1}{50-t} \cdot \frac{50-t}{3} = \frac{1}{3}\)

fj1 <- μ1/μτ; fj1[1]
## [1] 0.3333333

\(f_{J|T}(2|t) = \frac{\mu_{x}^{(2)}(t)}{\mu_{x}^{(\tau)}(t)} = \frac{2}{50-t} \cdot \frac{50-t}{3} = \frac{2}{3}\)

fj2 <- μ2/μτ; fj2[1]
## [1] 0.6666667

e) Prepare uma tabela e desenhe o gráfico de \({_t}p_{x}^{(\tau)}\) para \(t = 1, 2, . . . , 10\) e para \(x = 50\).

t = seq(1, 10, 1)
idades <- seq(51, 60, 1)
idades
##  [1] 51 52 53 54 55 56 57 58 59 60
dfp10τ <- data.frame(idades, pτ[2:11])
dfp10τ
##    idades pτ.2.11.
## 1      51 0.941192
## 2      52 0.884736
## 3      53 0.830584
## 4      54 0.778688
## 5      55 0.729000
## 6      56 0.681472
## 7      57 0.636056
## 8      58 0.592704
## 9      59 0.551368
## 10     60 0.512000
plot(idades, pτ[2:11], "o", pch = 20, cex=0.7, main="Gráfico",xlab="x+t", ylab="10p50(τ)")

f) Prepare uma tabela e desenhe o gráfico de \({_t}q_{x}^{(1)}\) e \({_t}q_{x}^{(2)}\) para \(t = 1, 2, . . . , 10\) e para \(x = 50\).

f1 <- function(t)(50-t)^2/(50^3)
f2 <- function(t)2*(50-t)^2/(50^3)

q11 <- rep(0, 10)
q12 <- rep(0, 10)
for (i in t) {
q11[i] <- integrate(f1, 0, i)$value
q12[i] <- integrate(f2, 0, i)$value
}
dfq10 <- data.frame(t, q11, q12); dfq10
##     t        q11        q12
## 1   1 0.01960267 0.03920533
## 2   2 0.03842133 0.07684267
## 3   3 0.05647200 0.11294400
## 4   4 0.07377067 0.14754133
## 5   5 0.09033333 0.18066667
## 6   6 0.10617600 0.21235200
## 7   7 0.12131467 0.24262933
## 8   8 0.13576533 0.27153067
## 9   9 0.14954400 0.29908800
## 10 10 0.16266667 0.32533333
plot(t, q11, type="o", pch=20, cex=0.7, col="blue", main=" ", xlab="t", ylab=bquote(q[x]), ylim=c(0, 0.35))
lines(t, q12, type="o", pch=20, cex=0.7, col="red")
legend("topleft", legend=c(bquote(q[50]^{(1)}), bquote(q[50]^{(2)})), col=c("blue", "red"), pch=20)

g) Encontre a distribuição de J(x) para x = 50.

\(f_{J} = {_t}q_{x}^{(j)}\)
\({_t}q_{50}^{(j)} = \int_{0}^{50}f_{T, J}(t, j) dt\)
\(_{50}q_{50}^{(1)} = \int_{0}^{50}\frac{(50 - t)^2}{50^3} dt\)

funcaoq501 <- function(t) (50 - t)^2 / (50^3)
q501 <- integrate(funcaoq501, 0, 50)$value; q501
## [1] 0.3333333

\(_{50}q_{50}^{(2)} = \int_{0}^{50} \frac{2 \cdot (50 - t)^2}{50^3} dt\)

funcaoq502 <- function(t) 2 * (50 - t)^2 / (50^3)
q502 <- integrate(funcaoq502, 0, 50)$value; q502
## [1] 0.6666667

h) Encontre a distribuição condicional de J(x) dado T(x) = 10 para x = 50.

\(f_{J|T}(1|10) = \frac{\mu_{50}^{(1)}(10)}{\mu_{50}^{(\tau)}(10)} = \frac{1}{50-10} \cdot \frac{50-10}{3} = \frac{1}{3}\)
\(f_{J|T}(2|10) = \frac{\mu_{50}^{(2)}(10)}{\mu_{50}^{(\tau)}(10)} = \frac{2}{50-10} \cdot \frac{50-10}{3} = \frac{2}{3}\)

2º - Elabore uma tábua de serviço para um grupo de 100 funcionários da empresa X, que participam do seu fundo de pensão. Considere que todos tenham a idade de 30 anos, o plano tenha taxa de juros de 6% a.a. e as seguintes taxas decrementais como hipóteses demográficas do plano:

I. Mortalidade a tabela: AT-2000 Male
II. Entrada em invalidez: ALVARO VINDAS;
III. Rotatividade: função monótona não crescente por idade no valor de 2% até 50 anos;
IV. Aposentadoria, segundo as probabilidades assumidas pela entidade fechada de previdência complementar, no período de 1998 a 2003:

Idade Probabilidade
55 0.10000
56 0.03000
57 0.03000
58 0.03000
59 0.03000
60+ 1.00000 
#importando tabuas de invalidez e mortalidade
tabua_invalidez <- read_excel("ALVARO_VINDAS.xlsx")
tabua_vida <- read_excel('AT-2000-male.xlsx')

#unindo as tabuas em um unico df
tabua_servico <- merge(tabua_vida, tabua_invalidez, by = "x")

#acrescentando a coluna de rotatividade
tabua_servico$rotatividade <- ifelse(tabua_servico$x >= 30 & tabua_servico$x <= 50, 0.02, 0)

#acrescentando a coluna de aposentadoria
tabua_servico$aposentadoria <- ifelse(tabua_servico$x == 55, 0.1, ifelse(tabua_servico$x >= 56 & tabua_servico$x <= 59, 0.03,
ifelse(tabua_servico$x >= 60, 1.0, 0)))

#filtrando a tabua de servico entre 30 e 60 anos
tabua_servico <- tabua_servico[tabua_servico$x >= 30 & tabua_servico$x <= 60, ]

# q_d taxa de mortalidade
tabua_servico$q_d <- tabua_servico$`mortalidade` * (1 - 1/2 * (tabua_servico$`invalidez` + tabua_servico$rotatividade + tabua_servico$aposentadoria) + 1/3 * (tabua_servico$`invalidez` * tabua_servico$rotatividade + tabua_servico$`invalidez` * tabua_servico$aposentadoria + tabua_servico$rotatividade * tabua_servico$aposentadoria) - 1/4 * tabua_servico$`invalidez` * tabua_servico$rotatividade * tabua_servico$aposentadoria)

# q_i taxa de invalidez
tabua_servico$q_i <- tabua_servico$`invalidez` * (1 - 1/2 * (tabua_servico$`mortalidade` + tabua_servico$rotatividade + tabua_servico$aposentadoria) + 1/3 * (tabua_servico$`mortalidade` * tabua_servico$rotatividade + tabua_servico$`mortalidade` * tabua_servico$aposentadoria + tabua_servico$rotatividade * tabua_servico$aposentadoria) - 1/4 * tabua_servico$`mortalidade` * tabua_servico$rotatividade * tabua_servico$aposentadoria)

# q_r taxa de rotatividade
tabua_servico$q_r <- tabua_servico$rotatividade * (1 - 1/2 * (tabua_servico$`mortalidade` + tabua_servico$`invalidez` + tabua_servico$aposentadoria) + 1/3 * (tabua_servico$`mortalidade` * tabua_servico$`invalidez` + tabua_servico$`mortalidade` * tabua_servico$aposentadoria + tabua_servico$`invalidez` * tabua_servico$aposentadoria) - 1/4 * tabua_servico$`mortalidade` * tabua_servico$`invalidez` * tabua_servico$aposentadoria)

# q_a taxa de aposentadoria
tabua_servico$q_a <- tabua_servico$aposentadoria * (1 - 1/2 * (tabua_servico$`mortalidade` + tabua_servico$`invalidez` + tabua_servico$rotatividade) + 1/3 * (tabua_servico$`mortalidade` * tabua_servico$`invalidez` + tabua_servico$`mortalidade` * tabua_servico$rotatividade + tabua_servico$`invalidez` * tabua_servico$rotatividade) - 1/4 * tabua_servico$`mortalidade` * tabua_servico$`invalidez` * tabua_servico$rotatividade)

# q_t taxa total
tabua_servico$q_t <- tabua_servico$q_d + tabua_servico$q_i + tabua_servico$q_r + tabua_servico$q_a
tabua_servico$p_t <- 1 - tabua_servico$q_t

# criando colunas lx_t e dx_t
for (i in 1:nrow(tabua_servico)) {
  if (i == 1) {
    tabua_servico$lx_t[i] <- 100
    tabua_servico$dx_t[i] <- tabua_servico$lx_t[i] * tabua_servico$q_t[i]
  } else {
    tabua_servico$lx_t[i] <- tabua_servico$lx_t[i-1] - tabua_servico$dx_t[i-1]
    tabua_servico$dx_t[i] <- tabua_servico$lx_t[i] * tabua_servico$q_t[i]
  }
}

# criando colunas de comutação Dx e Cx
v <- 1/(1+0.06)
tabua_servico$Dx <- tabua_servico$lx_t * v^tabua_servico$x
tabua_servico$Cx <- tabua_servico$dx_t * v^(tabua_servico$x+1)
for (i in 1:nrow(tabua_servico))
{
tabua_servico$Mx[i] = sum(tabua_servico$Cx[i:nrow(tabua_servico)])
}

tabua_servico <- round(tabua_servico, 6)
#tabua_servico
datatable(tabua_servico, options = list(pageLength = 31, scrollY = "300px", scrollX = TRUE))

a) Calcule a probabilidade de entrada em invalidez para a idade de 45 anos.

q_i_45 <- tabua_servico$q_i[tabua_servico$x == 45]
q_i_45
## [1] 0.001161

b) O número de ativos que provavelmente entrarão em invalidez aos 45 anos e irão sobreviver.

n_i_45 <- tabua_servico$lx_t[tabua_servico$x == 45] * tabua_servico$q_i[tabua_servico$x == 45] * (1 - tabua_servico$q_d[tabua_servico$x == 45]/2)
n_i_45
## [1] 0.08342029

c) Para um indivíduo de idade x = 48 anos, calcule o prêmio único de um seguro temporário de 10 anos (B = 1) que cobre a sua saída do fundo de pensão por todas as causas.

pup <- (tabua_servico$Mx[tabua_servico$x == 48] - tabua_servico$Mx[tabua_servico$x == 58]) / tabua_servico$Dx[tabua_servico$x == 48]

pup <- paste("RS", format(pup, decimal.mark = ",", big.mark = ".", nsmall = 2, digits = 6))

pup
## [1] "RS 0,17791"