Afleveringliv
Oskar Allerslev, Hugo Hendil, Emanuel Konrad
2024-12-11
Indlæsning af data
m1 <- read.csv("mort1.csv", sep = ";", dec = ",")
m1fun <- approxfun(m1$Alder, m1$mu, rule = 2)
mu02star <- function(t){ m1fun(t+30) }Opgave a)
I opgave a) og b) vil vi undertrykke notationen for reserverne og de tilhørende objekter.
Vi bemærker at der ikke er klumpbetalinger. Vi opskriver vores reward som \(R(t)=\Lambda^1(t) * B(t)+ \Delta(b(t))\).
\[ \begin{align*} R^\circ(t)=\begin{bmatrix} \left(-\pi \cdot 1_{\{t\leq 35 \}}+400K \cdot 1_{\{t > 35\}}\right) & 0 & 0 \\ 0 & 400K & 0 \\ 0 & 0 & 0 \end{bmatrix} \end{align*} \]
Vi har dertil også:
\[ \begin{align*} dB^\circ(t) &= \sum_{i \in E} \left( 1_{\{ X(t-)=i\}} b^\circ_i(t)dt\right)\\ &=-\pi \cdot 1_{\{t\leq 35 \cap X(t-)=1 \}}dt+400K \cdot 1_{\{t> 35 \cap X(t-)=1 \}}dt+400K \cdot 1_{\{ X(t-)=2 \}}dt \end{align*} \] Det skyldes vi ikke har klumpbetalinger. Den afhænger ikke af antagelserne fra førsteordensgrundlaget.
opgave b)
En vigtig antagelse for at man kan regne ækvivalenspræmien er at \(R(t)=R(t;\theta)\). Det er også gældende i vores tilfælde, da \(\theta\) bare er præmieraten, hvorfor vi kan benytte os af Newtons algoritme til at finde ækvivalenspræmien. Vi vælger at skrive vores udtryk på vektor-form. Bemærk \(t\in(0,80)\). Dertil er de prospektive reserver givet:
\[ \begin{align*} V^{th}(t;\pi) &= \int_{t}^{80}e^{-\int_{t}^u 0.01dx} P(0,u)R^{\circ *}(u)edu\\ &=\int_{t}^{80} e^{-\int_t^u 0.01 dx} \begin{bmatrix} p_{11}(0, u) & p_{12}(0, u) & p_{13}(0, u) \\ p_{21}(0, u) & p_{22}(0, u) & p_{23}(0, u) \\ p_{31}(0, u) & p_{32}(0, u) & p_{33}(0, u) \end{bmatrix} \begin{bmatrix} \left(-\pi \cdot 1_{\{t\leq 35 \}}+400K \cdot 1_{\{t > 35\}}\right) \\ 400K \\ 0 \end{bmatrix} du \\ &=\begin{bmatrix} \int_t^{80}e^{-\int_t^u 0.01 dx}p_{11}(0,u)\left(-\pi \cdot 1_{\{t\leq 35 \}}+400K \cdot 1_{\{t > 35\}}\right)+p_{12}(0,u)(400K )du\\ \int_t^{80}e^{-\int_t^u 0.01 dx}p_{21}(0,u)\left(-\pi \cdot 1_{\{t\leq 35 \}}+400K \cdot 1_{\{t > 35\}}\right)+p_{22}(0,u)(400K )du\\ 0 \end{bmatrix}. \end{align*} \] Vi har også \[ \begin{align*} V^{th}_1(0;0)=\int_0^{80}e^{-\int_0^u 0.01 dx}p_{11}(0,u)\left(400K \cdot 1_{\{t > 35\}}\right)+p_{12}(0,u)(400K )du, \end{align*} \] samt \[ \begin{align*} \frac{\partial}{\partial \pi}V^{th}_1(0;\pi)= \int_0^{80}e^{-\int_0^u 0.01 dx}p_{11}(0,u)\left(-\pi \cdot 1_{\{t\leq 35 \}}\right)du \end{align*} \] Hertil er ækvivalenspræmien givet ved: \[ \begin{align*} \mu = \frac{V^{th}_1(0;0)}{\frac{\partial}{\partial \pi}V^{th}_1(0;\pi)} = \frac{\int_0^{80}e^{-\int_0^u 0.01 dx}p_{11}(0,u)\left(400K \cdot 1_{\{t > 35\}}\right)+p_{12}(0,u)(400K )du }{ \int_0^{80}e^{-\int_0^u 0.01 dx}p_{11}(0,u)\left(-\pi \cdot 1_{\{t\leq 35 \}}\right)du} \end{align*} \]
Opgave c)
Vi benytter RK4 til at regne vores produktintegral. Nedenstående kode er (i stort omfang) kopieret fra den kode vi har fået udleveret. Den udgør vores beregning af produktintegralet.
prodint <- function(A, s, t, n) {
x0 <- s
y0 <- diag(nrow(A(s)))
h <- (t - s) / n
for (i in 1:n) {
s1 <- h * y0 %*% A(x0)
s2 <- h * (y0 + s1 / 2) %*% A(x0 + h / 2)
s3 <- h * (y0 + s2 / 2) %*% A(x0 + h / 2)
s4 <- h * (y0 + s3) %*% A(x0 + h)
y0 <- y0 + s1 / 6 + s2 / 3 + s3 / 3 + s4 / 6
x0 <- x0 + h
}
return(y0)
}
Lambda <- function(x) {
A <- matrix(0, 3, 3)
A[1, 2] <- (0.0004 + 10^(4.54 + 0.06*(x+30)-10))*ifelse(x <= 35, 1, 0)
A[1,3] <- mu02star(x)
A[2,1] <- 2.0058 * exp(-0.117*(x+30)) * ifelse(x <= 35,1,0)
A[2,3] <- A[1,3]*(1+ifelse(x <= 35,1,0))
row_sums <- rowSums(A)
diag(A) <- -row_sums
return(A)
}
R <- function(x, mu) {
if (x <= 35) {
return(diag(c(-mu , 400000, 0)))
} else {
return(diag(c(400000, 400000,0)))
}
}
dR <- function(x, mu) {
if (x <= 35) {
return(diag(c(-1, 0,0)))
} else {
return(diag(c(0,0,0)))
}
}
reserve <- function(t, TT, Lambda, R, mu, rint, N) {
dim <- nrow(Lambda(t))
A11 <- function(x) {
return(Lambda(x) - rint * diag(1, nrow = dim))
}
RM <- function(x) {
cbind(rbind(A11(x), matrix(0, dim, dim)), rbind(R(x, mu), Lambda(x)))
}
PRM <- prodint(RM, t, TT, N)
RES <- PRM[1:dim, (dim + 1):(2 * dim)]
return((RES %*% rep(1, dim)))
}
equiv_premium <- function(t, TT, Lambda, R, dR, mu, rint, numsteps) {
b <- reserve(t, TT, Lambda, R, 0.0, rint, numsteps)
a <- reserve(t, TT, Lambda, dR, 0.0, rint, numsteps)
return(-b / a)
}
equiv_premium(0, 80, Lambda, R, dR, 0.0, 0.01, 5000)[1]## [1] 209640.5Opgave d)
V_20_ <- reserve(20, 80, Lambda, R, 209640.5, 0.01, 1000)
V_35_ <- reserve(35, 80, Lambda, R, 209640.5, 0.01, 1000)
data.frame( "V(20)" = V_20_, "V(35)" = V_35_)Vi observerer vores resultater er i overenstemmelse med den standard fortolkning af prospektive reserver.
opgave e)
Vi benytter \[ \begin{align*} \frac{\partial}{\partial t} V_i^{\circ, *}(t)=(r^*(t)-\Lambda^*(t))V_i^{\circ, *}(t)-R^{\circ *}(t)e. \end{align*} \]
TT <- 80
eq <- 209640.5
rint <- 0.01
N <- 1000
x_vals <- seq(0, TT, length.out = 100)
y_vals <- sapply(x_vals, function(t) {
res_mat <- reserve(t, TT, Lambda, R, eq, rint, N)
v <- res_mat
c(V_1 = v[1], V_2 = v[2])
})
y_vals <- t(y_vals)
results <- data.frame(x = x_vals, V_1 = y_vals[,1], V_2 = y_vals[,2])
par(mfrow = c(1, 2))
plot(results$x, results$V_1, xlab = 't', ylab = 'Reserve', main = 'V_1 Reserve')
abline(0,0)
plot(results$x, results$V_2, xlab = 't', ylab = 'Reserve', main = 'V_2 Reserve')
abline(0,0)
Vi bemærker at reserverne for de to tilstande er som forventet. De er i overenstemmelse med den standard fortolkning af prospektive reserver for henholdsvis aktiv og invalid. ‘Knækket’ omkring \(t=70\) kan skyldes, vi ikke skelner mellem intensiteter for mænd og kvinder. Det kan også tænkes at være en måde vi inkorporere halerisikoen for de ældre.
diff_partial_thiele <- function(t,TT=80, eq =209640.5 , rint = 0.01, N=1000){
result <- (rint - Lambda(t)) %*% reserve(t, TT, Lambda, R, eq, rint, N) - R(t,eq) %*% rep(1,3)
c(V_1 = result[1], V_2 =result[2])
}
x_vals <- seq(0, 80, length.out = 100)
y_vals <- t(sapply(x_vals, diff_partial_thiele))
results <- data.frame(x = x_vals, V_1 = y_vals[, 1], V_2 = y_vals[, 2])
plot(results$x, results$V_1, type = "l", col = "blue", lwd = 2,
xlab = "x (t values)", ylab = "Values", main = "V_1 and V_2 vs x")
abline(a = 0, b = 0)
lines(results$x, results$V_2, col = "red", lwd = 2)
legend("topright", legend = c("Aktiv", "Invalid"), col = c("blue", "red"), lty = 1, lwd = 2)
Man kan igen her se at differentialligningerne er som forventet og i
overenstemmelse med de reserver vi plottede ovenfor denne.
Opgave f)
Vi finder her standardafvigelsen ved brug af \((2.29)\) (udtrykket i forelæsningsnoterne) til at finde \((2.30)\) (udtrykket i forlæsningsnoterne) hvilket konkret beregnes ved: \[ \begin{align*} \prod_s^t (I+F^{(k)}(x)dx)=H^{(k)}(s,t). \end{align*} \]
eq <- 209640.5
R <- function(x, mu) {
if (x <= 35) {
return(diag(c(-mu , 400000, 0)))
} else {
return(diag(c(400000, 400000,0)))
}
}
Fu <- function(Lambda, R){
I <- base::diag(x=1, nrow = 3, ncol = 3)
k=2
r = 0.01
A <- matrix(0,9,9)
A[1:3, 1:3] <- Lambda- k*r*I
A[4:6, 4:6] <- Lambda - (k-1)*r*I
A[7:9, 7:9] <- Lambda - (k-2)*r*I
A[1:3, 4:6] <- R
A[4:6,7:9] <- R
return(A)
}
ABE <- function(u){
Fu(Lambda(u),R(u,eq))
}
resultat <- prodint(ABE, 0, 80, 10000)
sqrt((2*resultat[1:3,7:9]-resultat[4:6,7:9]^2) %*% rep(1,3))[1]## [1] 2952402Opgave g)
Først definerer vi vores nye rentefunktion vha. interpolation (som også gjort tidligere). Derudover konstruerer vi også vores \(\Lambda(t)\).
m2 <- read.csv("ftmort.csv", sep = ";", dec = ",")
m2fun <- approxfun(m2$Alder, m2$mu, rule = 2)
mu202star <- function(t){ m2fun(t+30) }
rente <- read.csv("rente.csv", sep = ";", dec = ",")
rentefun <- approxfun(rente$t, rente$rt, rule = 2)
rentefunk <- function(t){rentefun(t)}
Lambda1 <- function(x) {
A <- matrix(0, 3, 3)
A[1, 2] <- (10^(5.662015+0.033462*(x+30)-10))*ifelse(x <= 35, 1, 0)
A[1,3] <- mu202star(x)
A[2,1] <- 4.0116*exp(-0.117*(x+30))*ifelse(x<=35,1,0)
A[2,3] <- (0.010339+10^(5.070927+0.05049*(x+30)-10))*ifelse(x <= 35,1 ,0 )+ A[1,3]*ifelse(x>35,1,0)
row_sums <- rowSums(A)
diag(A) <- -row_sums
return(A)
}Vi regner nu risikosummerne. Jf. vi stadig ikke har nogle klumpbetalinger benyttes: \[ \begin{align*} RS^{*, \circ}_{ik}(t)=V_k^{*, \circ}(t)-V_i^{*, \circ}(t). \end{align*} \]
TT <- 80
eq <- 209640.5
rint <- 0.01
N <- 1000
time_points <- seq(0, TT, length.out = 100)
differences <- numeric(length(time_points))
for (i in seq_along(time_points)) {
t <- time_points[i]
fo_thiele <- reserve(t, TT, Lambda, R, eq, rint, N) %*% rep(1, 3)
differences[i] <- fo_thiele[2] - fo_thiele[1]
}
plot(time_points, differences, type = "l", col = "blue", lwd = 2,
xlab = "Time (t)", ylab = "Difference (V1 - V2)",
main = "Difference between fo_thiele[1] and fo_thiele[2] over Time")
abline(a = 0, b = 0)
Vi regner nu det samme for tredjeordensgrundlaget:
reserveR <- function(t, TT, Lambda, R, mu, r, N) {
dim <- nrow(Lambda(t))
A11 <- function(x) {
return(Lambda(x) - r(x) * diag(1, nrow = dim))
}
RM <- function(x) {
cbind(rbind(A11(x), matrix(0, dim, dim)), rbind(R(x, mu), Lambda(x)))
}
PRM <- prodint(RM, t, TT, N)
RES <- PRM[1:dim, (dim + 1):(2 * dim)]
return((RES %*% rep(1, dim)))
}
time_points_1 <- seq(0, 80, length.out = 100)
differences1 <- numeric(length(time_points))
for (i in seq_along(time_points)) {
t <- time_points[i]
fo_thiele <- reserveR(t,80,Lambda1,R,eq,rentefunk,1000)
differences1[i] <- fo_thiele[2] - fo_thiele[1]
}
plot(time_points, differences, type = "l", col = "blue", lwd = 2,
xlab = "Time (t)", ylab = "Difference (V1 - V2)",
main = "Difference between fo_thiele[1] and fo_thiele[2] over Time")
lines(time_points_1, differences1, col = "red", lwd = 2)
Derfra har vi \(RS_{1,2} \leq
RS_{1,2}^{*,\circ}\). De resterende tilfælde tolker vi som
unødige at regne.
Opgave h)
Vi vil nu finde transitionssandsynlighederne for vores markovkæde. Vi benytter os af vores produktintegral til at finde disse. Vi benytter os af følgende algoritme.
start_time <- 0
end_times <- 0:80
num_steps <- 1000
devtools::install_github("OskarAllerslev/lifepack")## WARNING: Rtools is required to build R packages, but is not currently installed.
##
## Please download and install Rtools 4.4 from https://cran.r-project.org/bin/windows/Rtools/.## Using GitHub PAT from the git credential store.## Skipping install of 'lifepack' from a github remote, the SHA1 (6fe00c46) has not changed since last install.
## Use `force = TRUE` to force installationlibrary(lifepack)##
## Vedhæfter pakke: 'lifepack'## De følgende objekter er maskerede _af_ '.GlobalEnv':
##
## equiv_premium, prodint, reserveP11 <- numeric(length(end_times))
P12 <- numeric(length(end_times))
P13 <- numeric(length(end_times))
for (i in seq_along(end_times)) {
end_time <- end_times[i]
prod_matrix <- lifepack::prodint(Lambda, start_time, end_time, num_steps)
P11[i] <- prod_matrix[1, 1]
P12[i] <- prod_matrix[1, 2]
P13[i] <- prod_matrix[1, 3]
}
plot(end_times, P11, type = "l", col = "blue", lwd = 2,
xlab = "Time", ylab = "Transition Probability",
main = "Transition Probabilities Over Time")
lines(end_times, P12, col = "red", lwd = 2)
lines(end_times, P13, col = "green", lwd = 2)
legend("bottomleft", legend = c("P(1 -> 1)", "P(1 -> 2)", "P(1 -> 3)"),
col = c("blue", "red", "green"), lty = 1, lwd = 2)
list(
data.frame(transition = "s = 25",
p11 = prodint(Lambda1, 0 ,25, 1000)[1,1],
p12 = prodint(Lambda1, 0 ,25, 1000)[1,2],
p13 = prodint(Lambda1, 0 ,25, 1000)[1,3]),
data.frame(transition = "s = 35",
p11 = prodint(Lambda1, 0 ,35, 1000)[1,1],
p12 = prodint(Lambda1, 0 ,35, 1000)[1,2],
p13 = prodint(Lambda1, 0 ,35, 1000)[1,3])
)## [[1]]
## transition p11 p12 p13
## 1 s = 25 0.9519117 0.02685201 0.02123634
##
## [[2]]
## transition p11 p12 p13
## 1 s = 35 0.8736452 0.05891422 0.06744058Vores resultater virker fornuftige.
Opgave i)
Vi regner det forventede tilstandsvise cash flow som givet ift.: \[ \begin{align*} a(s,t)=P(s,t)R(t) e. \end{align*} \]
Dette giver en vektor, hvoraf vi kan hive resultatet for tilstand \(i = 1,2,3\) ud.
eq <- 209640.5
start <- 0
end_s <- 0:80
steps <- 1000
e_vec <- c(1, 1, 1)
a0_0 <- numeric(length(end_s))
a0_1 <- numeric(length(end_s))
for (i in seq_along(end_s)) {
s <- end_s[i]
P_mat <- prodint(Lambda, start, s, steps)
R_mat <- R(s, eq)
a_vec <- P_mat %*% R_mat %*% e_vec
a0_0[i] <- a_vec[1]
a0_1[i] <- a_vec[2]
}
plot(end_s, a0_0, type = "l", col = "blue", lwd = 2,
xlab = "Time s", ylab = "Expected cash flow")
lines(end_s, a0_1, col = "red", lwd = 2)
abline(0,0)
sel_s <- c(20, 35)
df_values <- data.frame(
s = sel_s,
a0_0 = a0_0[end_s %in% sel_s],
a0_1 = a0_1[end_s %in% sel_s]
)
df_valuesOpgave j)
Vi skal regne markedsværdien for kontrakten, hvilket er reserven udregnet på baggrund af markedsværdigrundlaget. \[ \begin{align*} V_0^\circ(0)=\int_0^{80}e^{\int_t^u r(v)dv}a_0^\circ(t,u)du. \end{align*} \]
eq <- 209640.5
t <- 0
reserveR(t,80, Lambda1, R, eq, rentefun, 1000)[1]## [1] -2160684opgave k)
\(\textbf{Antag}\): \(c(t) - \delta(t) > 0 \hspace{0.5cm} \forall t \in \mathbb{R}_+\).
Vi ved at \(V^*(t) \geq 0 \hspace{0.5cm} \forall t \in \mathbb{R}_+\), da vi benytter ækvivalenspræmien regnet ud fra førsteordensgrundlaget. Fra tidligere inspektion af Risikosummerne har vi at
\[ \begin{align*} RS_{01}(t) \geq & 0\\ RS_{02}(t) \leq & 0\\ RS_{10}(t) \leq & 0 \\ RS_{12} \leq & 0\\ \end{align*} \] Dertil kan vi ved benyttelse af theorem 6.5 skrive: \[ \begin{align*} c(t) = (r(t)-r^\delta(t))V^*(t) + \sum_{j \in E \land j \neq X(t, \omega)} (\mu_{X(t,\omega),j}(t)- \mu_{X(t,\omega),j}(t))RS^*_{X(t,\omega),j} > \delta(t) \end{align*} \] Hvorved vi benytter vores antagelse om at \(c(t) - \delta(t) > 0\) og derved udleder de nødvendige krav for \(r^\delta(t)\) og \(\alpha_{i,j}\). For at ovenstående skal gælde har vi
\[ \begin{align*} r(t) & \geq r^\delta(t)\\ \mu_{01}^\delta(t) & \geq \mu_{01}(t)\\ \mu_{02}^\delta(t) & \leq \mu_{02}(t)\\ \mu_{10}^\delta(t) & \leq \mu_{10}(t)\\ \mu_{12}^\delta(t) & \leq \mu_{12}(t) \end{align*} \]
Vi kan nu benytte os af vores \(\mu_{01}^\delta(t)\) og \(\mu_{02}^\delta(t)\) til at finde vores \(\alpha_{i,j}\).
\[ \begin{align*} \alpha_{01} = \frac{\mu_{01}^\delta(t)}{\mu_{01}(t)} \geq 1\\ \alpha_{02} = \frac{\mu_{02}^\delta(t)}{\mu_{02}(t)} \leq 1\\ \alpha_{10} = \frac{\mu_{10}^\delta(t)}{\mu_{10}(t)} \leq 1\\ \alpha_{12} = \frac{\mu_{12}^\delta(t)}{\mu_{12}(t)} \leq 1 \end{align*} \]
hvorved vi har fundet de nødvendige krav.