Chapter 3

###3.1 Market portfolio and CAPM

###3.1.1 Mean-variance portfolio

rA<- c(-0.2,0.1,0.3,0.5)
rB<- c(0.05,0.2,-0.12,0.09)
prob <- rep(0.25,4)

expA <-sum(rA * prob)

expA

expB <-sum(rB *prob)

expB

varA <-sum(rA^2*prob)-expA^2
varA

varB <-sum(rB^2*prob)-expB^2
varB

covAB <- sum(rA*rB*prob)-expA*expB

covAB

#[Link] expected return and the variance of portfolio as:

wA<- 0.6
wB<- 1-wA

exp_por <- wA*expA+wB*expB

exp_por

var_por <- varA*wA^2+varB*wB^2+2*covAB*wA*wB

var_por

#[Link] standard deviation（小）:

sqrt(var_por)

#[Link] weighted average of the standard deviations of RA and RB（大）:

wA*sqrt(varA)+wB*sqrt(varB)

#which is greater than the standard deviation of the portfolio,this is known as the diversification
effect（多样化效应）
#另一种方法求解 exp_pro and sd_pro
w <- c(wA, 1 - wA)
exp_por <- w %*% c(expA, expB)
exp_por
sigma <- matrix(c(varA, covAB, covAB, varB), ncol = 2)
sd_pro <- sqrt(t(w) %*% sigma %*% w)
sd_pro

#5 opportunity set
wa <- seq(0, 1, by = 0.01)
w <- cbind(wa, 1 - wa)
exp_p <- w %*% c(exp_A, exp_B)
sigma <- matrix(c(var_A, covAB, covAB, var_B), ncol = 2)
sd_portfolio <- function(sigma, w) {
sg <- rep(0, nrow(w))
for (i in 1:nrow(w)) {
sg[i] <- sqrt(t(w[i, ]) %*% sigma %*% (w[i, ]))
}
sg
}
sd_p <- sd_portfolio(sigma, w)

plot(sd_p, exp_p,
type = "l",
main = "Opportunity set",
xlab = "SD(Rp)",
ylab = "E(Rp)"
)
#get the plot

#6 minimum variance portfolio（MV）

在 opportunity set 上，有一个点使 SD(Rp)最小值，称为最小方差投资组合（MV）
wmin <- w[match(min(sd_p), sd_p), ]
wmin
#最小方差投资组合对应的期望收益 exp_mv 和标准差 sd_mv：
exp_mv <- wmin %*% c(expA, expB)
exp_mv

sd_mv <- sd_portfolio(sigma, matrix(wmin, 1))

sd_mv

假设我们想要通过混合资产来创建一个投资组合，使得最终投资组合的标准差大于可以实现的最小标准差
我们有两个选择在一个相同的标准差上，通常选择 E【Rp】较高的投资组合
因此，在点 MV 上方的机会集被称为有效边界

#7 efficient frontier 有效边界

exp_ef <- exp_p[exp_p > [Link](exp_mv)]
sd_ef <- sd_p[exp_p > [Link](exp_mv)]
plot(sd_p, exp_p,
type = "l",
main = "Opportunity set",
xlab = "SD(Rp)",
ylab = "E(Rp)"
)
points(sd_mv, exp_mv, col = "red", pch = 19) # MV
lines(sd_ef, exp_ef, col = "blue", lwd = 2) # Efficient frontier

#8 机会集的形状是如何根据风险资产的相关性变化的？
Question:or that purpose, let us consider two stocks, A and B, with the following characteristic:
E [ R A ] = 10%, E [ R B ] = 15%, SD (R A ) = 10%, and SD (R B ) = 12%. We also assume five di fferent
possibilities for the correlation between A and B: -1, -0.5, 0, 0.5, and 1

wa <- seq(0, 1, by = 0.01)

w <- cbind(wa, 1 - wa)
exp_p <- w %*% c(0.1, 0.15)
sd_A <- 0.1
sd_B <- 0.12
sigma1 <- matrix(c(sd_A^2, -1 * sd_A * sd_B, -1 * sd_A * sd_B, sd_B^2), 2)
sigma2 <- matrix(c(sd_A^2, -0.5 * sd_A * sd_B, -0.5 * sd_A * sd_B, sd_B^2), 2)
sigma3 <- matrix(c(sd_A^2, 0 * sd_A * sd_B, 0 * sd_A * sd_B, sd_B^2), 2)
sigma4 <- matrix(c(sd_A^2, 0.5 * sd_A * sd_B, 0.5 * sd_A * sd_B, sd_B^2), 2)
sigma5 <- matrix(c(sd_A^2, 1 * sd_A * sd_B, 1 * sd_A * sd_B, sd_B^2), 2)
sd_p1 <- sd_portfolio(sigma1, w)
sd_p2 <- sd_portfolio(sigma2, w)
sd_p3 <- sd_portfolio(sigma3, w)
sd_p4 <- sd_portfolio(sigma4, w)
sd_p5 <- sd_portfolio(sigma5, w)
plot(sd_p1, exp_p,
type = "l",
main = "Opportunity set",
xlab = "SD(Rp)",
ylab = "E(Rp)"
)
lines(sd_p2, exp_p, col = "red")
lines(sd_p3, exp_p, col = "blue")
lines(sd_p4, exp_p, col = "green")
lines(sd_p5, exp_p, col = "orange")
legend("topleft",
leg = paste("rho =", c(-1, -0.5, 0, 0.5, 1)),
lty = 1,
col = c("black","red","blue","green","orange")
)
Example 3.1. Consider three stocks, A, B, and C, with expected returns E [ R A ] = 0 .105, E [ R B ] =
0.18, and E[R C ] = 0.02, respectively. The covariance matrix is the following:
Then, the following code plots the opportunity set available to any investor

exp_r <- c(0.105, 0.18, 0.02)

sigma <- matrix(
c(
0.15^2, -0.012, 0.002,
-0.012, 0.12^2, -0.002,
0.002, -0.002, 0.2^2
),
ncol = 3, byrow = T
)
sigma
# Generates a grid of values for w1 and w2
w_grid <- [Link](
wa = seq(0, 1, [Link] = 100),
wb = seq(0, 1, [Link] = 100)
)
w_grid <- [Link](w_grid)
w <- cbind(w_grid, 1 - w_grid[, 1] - w_grid[, 2])
# Expected return
rp <- [Link](w %*% exp_r)
sd_p <- sd_portfolio(sigma, w)
plot(sd_p, rp,
pch = ".",
main = "Opportunity set",
xlab = "SD(Rp)",
ylab = "E(Rp)"
)

###3.1.2 Capital Asset Pricing Model (CAPM)

R = wpRp + wfRf = wpRp + (1 − wp)Rf .

# opportunity set of the (complete) portfolio

rf <- 0.12 # Risk-free rate
slope1 <- (rp[50] - rf) / sd_p[50]
slope2 <- (rp[2500] - rf) / sd_p[2500]
slope3 <- (rp[7500] - rf) / sd_p[7500]
x <- seq(0, 1, by = 0.01)
plot(sd_p, rp,
pch = ".",
xlab = "SD(R)",
ylab = "E(R)",
xlim = c(0, 0.23)
)
lines(x, slope1 * x + rf)
lines(x, slope2 * x + rf)
lines(x, slope3 * x + rf)
points(sd_p[50], rp[50], col = "orange", pch = 19)
points(sd_p[2500], rp[2500], col = "green", pch = 19)
points(sd_p[7500], rp[7500], col = "blue", pch = 19)

#the capital market line(CML) and market portfolio

rf <- 0.12 # Risk-free rate

amax <- max((rp - rf) / sd_p) # Slope of the CML
x <- seq(0, 1, by = 0.01)
plot(sd_p, exp_p,
pch = ".",
xlab = "SD(Rp)",
ylab = "E(Rp)",
xlim = c(0, max(sd_p))
)
lines(x, amax * x + rf)
points(sd_p[orp], rp[orp], col = "red", pch = 19)

#the weights of the market portfolio

orp <- match(amax, (rp - rf) / sd_p) # Optimal risky portfolio
w[orp, ]

unscaled_w <- solve(sigma, exp_r - rf)

scaled_w <- unscaled_w / sum(unscaled_w)
scaled_w

#assume that the following vectors represent the historical returns of the market portfolio (r_m) and
a risky asset (r_ra)
#假设以下向量代表市场投资组合（r_m）和风险资产（r_ra）的历史回报
r_m <- c(
0.000341, 0.068962, 0.046964, 0.006922,
-0.029561, 0.028035, -0.027218, -0.161576,
0.060479, 0.071397, 0.058284
)
r_ra <- c(
-0.050484, 0.194222, -0.021584, -0.009475,
0.070406, -0.074354, 0.078366, -0.193192,
-0.059271, 0.148521, 0.11004
)

#We now plot a graph of (r_m, r_ra) and draw the regression line that best fits the data
#绘制一个（r_m，r_ra）的图，并绘制出最适合数据的回归线

dat_beta <- [Link](r_m = r_m, r_ra = r_ra)

reg_line <- lm(r_ra ~ r_m, data = dat_beta)
plot(dat_beta)
abline(reg_line, col = "blue")

#The coefficients of this regression line 回归线的系数：

reg_line

#回归的斜率就是风险资产的贝塔系数。
#对 beta 的自然解释是，风险资产的回报率比市场的回报率放大了 1.047370 倍
#the returns of the risky asset are magnified 1.047370 times over those of the market.

###3.2 The binomial model

###3.2.1 One period binomial model

Example3.2
Consider a stock with a current price of S0 = 50. The stock’s price is expected to increase to 60 or
decrease to 40 during the next year. The risk-free interest rate is 5% compounded annually.
Compute the price of a 1-year European call option with a strike price of 55.

s0 <- 50
su <- 60
sd <- 40
rf <- 0.05
strike <- 55
u <- su / s0
d <- sd / s0
# Note that d <= (1 + R) <= u is satisfied
u
d
# Risk-neutral probabilities
qu <- (1 + rf - d) / (u - d)
qu
qd <- 1-qu
qd
# Value of the contingent claim at t = 1
phiu <- max(su - strike, 0)
phiu
phid <- max(sd - strike, 0)
phid
# Derivative price
price <- (phiu * qu + phid * qd) / (1 + rf)
price

###3.2.2 Multiperiod binomial model

Example 3.3.
Consider a stock with a current price of S 0 = 50. The stock’s price is expected to increase by 10% or
decrease by 8% during the next two six-month periods. The risk-free interest rate is 5%
compounded annually. Compute the price of a 1-year European put option with a strike price of
55.

s0 <- 50
u <- 1.1
d <- 0.92
strike <- 55
# Risk-free rate for the period of six months
rf <- (1 + 0.05)^(6 / 12) - 1
# Number of periods
n <- 2
# Future stock prices
k <- 0:n
s <- u^k * d^(n - k) * s0
s
# Value of the contingent claim at maturity
phi <- pmax(strike - s, 0)
phi
# Risk-neutral probability qu
qu <- (1 + rf - d) / (u - d)
qu
# Price of the option
vu <- (qu * phi[3] + (1 - qu) * phi[2]) / (1 + rf)
vu
vd <- (qu * phi[2] + (1 - qu) * phi[1]) / (1 + rf)
vd
v0 <- (qu * vu + (1 - qu) * vd) / (1 + rf)
v0

In other words, the binomial algorithm consists of the following steps:

1. Generate the tree of stock prices.
2. Calculate the value of the claim at each final node.
3. Calculate the claim values sequentially at each preceding node.
换句话说，二项式算法包括以下步骤：
1.生成股票价格的树。
2.计算每个最终节点上的索赔值。
3.在前面的每个节点上依次计算索赔值。

# To generate the tree of stock prices 生成股票价格的树

build_tree <- function(s0, u, d, n) {
tree <- matrix(0, nrow = n + 1, ncol = n + 1)
for (t in 1:(n + 1)) {
k <- 0:(t - 1)
tree[t, 1:t] <- u^k * d^(t - 1 - k) * s0
}
tree
}
tree <- build_tree(50, 1.1, 0.92, 2)
tree

#evaluate the claim

pmax(strike - tree[nrow(tree), ], 0)

#recursive computations
value_bin_mod <- function(qu, rf, tree, strike) {
val_tree <- matrix(0, nrow = nrow(tree), ncol = ncol(tree))
val_tree[nrow(tree), ] <- pmax(strike - tree[nrow(tree), ], 0) # European put
for (t in (nrow(tree) - 1):1) {
for (k in 1:t) {
val_tree[t, k] <- ((1 - qu) * val_tree[t + 1, k]
+ qu * val_tree[t + 1, k + 1]) / (1 + rf)
}
 }
val_tree
}
opt_price <- value_bin_mod(qu, rf, tree, strike)
opt_price

opt_price[1, 1]

#另一种解决方案：
pi0 <- sum(choose(n, k) * qu^k * (1 - qu)^(n - k) * phi) / (1 + rf)^n
pi0

R packages for the binomial model

library(derivmkts)
[Link](derivmkts)

#To compute prices of European (and American) call and put options 计算欧洲（美国）的看涨/看跌期权价格
we use the function binomopt()

binomopt(
s = 50, # Initial stock price
k = 55, # Strike price
r = log(1 + 0.05), # Continuously-compounded risk-free rate
tt = 1, # Time to maturity
d = 0, # Dividends, in our case, we do not work with dividends, hence 0.
nstep = 2, # Number of periods
american = FALSE, # To indicate European
putopt = TRUE, # To indicate a Put option
specifyupdn = TRUE, # Tells the function to use u and d
up = 1.1, # Value of u
dn = 0.92 # Value of d
)

# plots the development of the stock price and shows graphically the probability of being at each
node
#绘制股票价格的发展，并以图形方式显示在每个节点的概率
binomplot(
s = 50, k = 55, r = log(1 + 0.05), tt = 1, d = 0, nstep = 2,
american = FALSE, putopt = TRUE, specifyupdn = TRUE, up = 1.1, dn = 0.92,
v = 0, # A value of the volatility must be provided, although not used
plotarrows = TRUE, # Plots arrows that connect the nodes of the tree
plotvalues = TRUE # Plots the values of the stock prices at each node
)
#绿色和红色表示该选项是否在那里得到最佳执行（如果是则为绿色，如果不是则为红色）
#Example 3.4
The price of a stock is currently 40. Over each of the following three 5-month periods, it
is expected to go up by 10% or down by 5%. The risk-free interest rate is 6% per annum with
continuous compounding.
1. Use a three-step binomial model to find the value of a European style derivative that pays off
X = [max(K − ST , 0)]3, where ST is the stock price in 15 months and K = 44.
2. Use a three-step binomial model to find the value of a European-style derivative that pays off
X = min(ST^2, K), where ST is the stock price in 15 months and K = 2000.

#1.
s0 <- 40
kk <- 44
u <- 1.1
d <- 0.95
rf <- exp(0.06 * 5 / 12) - 1
tree <- build_tree(s0, u, d, 3)
tree

qu <- (1 + rf - d) / (u - d)
qu

value_bin_mod <- function(qu, rf, tree, kk) {

val_tree <- matrix(0, nrow = nrow(tree), ncol = ncol(tree))
val_tree[nrow(tree), ] <- (pmax(kk - tree[nrow(tree), ], 0))^3 # Given claim
for (t in (nrow(tree) - 1):1) {
for (k in 1:t) {
val_tree[t, k] <- ((1 - qu) * val_tree[t + 1, k]
+ qu * val_tree[t + 1, k + 1]) / (1 + rf)
}
}
val_tree
}
opt_price <- value_bin_mod(qu, rf, tree, kk)
opt_price

#Thus,the price of the option is:

opt_price[1, 1]

#2.
n <-3
k <- 0:n
# Prices after 15 months
s <- u^k * d^(n - k) * s0
s

# Value of the contingent claim at maturity

phi <- pmin(s^2, 2000)
phi

pi0 <- sum(choose(n, k) * qu^k * (1 - qu)^(n - k) * phi) / (1 + rf)^n

pi0

###3.3The Black and Scholes model

###3.2.1 Preliminars: Brownian motion

#模拟布朗运动的轨迹
[Link](1)
delta <- 0.001 # Increment
t <- seq(0, 1, by = delta) # Time interval
w <- rnorm(n = length(t) - 1, sd = sqrt(delta)) # iid normal distributed r.v.s
w <- c(0, cumsum(w)) # Cumulative sum - 0 is the initial value
plot(t, w,
type = "l",
xlab = "Time",
ylab = "W(t)",
main = "Simulated trajectory of a standard Brownian motion",
)

#生成多个轨迹
nsim <- 100 # Number of simulated trajectories
w <- matrix(
rnorm(n = nsim * (length(t) - 1), sd = sqrt(delta)),
nsim, # Each row is one simulated trajectory
length(t) - 1
)
w <- cbind(rep(0, nsim), t(apply(w, 1, cumsum))) # matrix with trajectories
plot(t, w[1, ], # Plots the first trajectory
type = "l",
ylim = c(-2.5, 2.5),
xlab = "Time",
ylab = "W(t)",
main = "Simulation of standard Brownian motion"
)
apply(w[-1, ], 1, function(x, t) lines(t, x), t = t) # Plots the remaining trajectories
#标准布朗运动可以推广到算术布朗运动，它是前者的尺度和移动。
#X(t)是一个算术布朗运动，X(t) = µt + σW(t)
# straightforward ways to simulate the arithmetic Brownian motion.

mu <- 2
sigma2 <- 0.2
nsim <- 100
x <- matrix(
rnorm(n = nsim * (length(t) - 1), mean = mu * delta, sd = sqrt(delta * sigma2)),
nsim,
length(t) - 1
)
x <- cbind(rep(0, nsim), t(apply(x, 1, cumsum)))
plot(t, x[1, ],
type = "l",
ylim = c(-1, 3),
xlab = "Time",
ylab = "X(t)",
main = "Simulation of arithmetic Brownian motion"
)
apply(x[-1, ], 1, function(x, t) lines(t, x), t = t)

#Note that the arithmetic Brownian motion can take negative values 算术布朗运动可以取负值
#S(t) = S(0) exp (X(t)) = S(0) exp (µt + σW(t))

s0 <- 2
s <- s0 * exp(x)
plot(t, s[1, ],
type = "l",
ylim = c(0, 40),
xlab = "Time",
ylab = "S(t)",
main = "Simulation of geometric Brownian motion"
)
apply(s[-1, ], 1, function(x, t) lines(t, x), t = t)

#3.3.2 The Black and Scholes formula

Example 3.5.
Consider a European call option over a stock with a current price of S(0) = 50 and
volatility of 0.25. Moreover, the risk-free interest rate with continuous compounding is 6% per
annum, the strike price is 45, and the option’s time to maturity is 6 months. Find the price of the
option.

s0 <- 50
sigma <- 0.25
rf <- 0.06
maturity <- 6 / 12
strike <- 45
d <- 0 # No dividends
bscall(s0, strike, sigma, rf, maturity, d)

# using the function simprice() to approximate an option’s price

[Link](1)
st <- simprice(s0, sigma, rf, maturity, d, trials = 10000, periods = 1)
exp(-rf * maturity) * mean(pmax(st[st$period == 1, ]$price - strike, 0))

#3.3.3 Greeks

#the Greeks for the option in Exercise 3.5

greeks(bscall(s0, strike, sigma, rf, maturity, d = 0))

#greeks() accepts vector inputs, which allows us to visualize the Greeks

s <- seq(.5, 80, by = .5)
call_greeks <- greeks(bscall(s, strike, sigma, rf, maturity, d = 0))
for (i in rownames(call_greeks)) {
plot(s, call_greeks[i, ], main = paste(i), ylab = i, type = "l", col = "red")
}

Chapter 4
###4.1 The collective risk model
###4.1.1 Discretization of claim amount distributions 索赔金额分配的自由裁量化

#Example 4.1.
Find the four discretizations of a Gamma(2, 1) distribution on (0, 10) with a step of 0.5 and plot
their cdf against the original cdf.

[Link]("actuar")
library("actuar")
# plot-upper
fx_upper <- discretize(pgamma(x, 2, 1),
method = "upper",
from = 0,
to = 10,
step = 0.5
)
x <- seq(0, 10 - 0.5, 0.5)
plot(stepfun(x, diffinv(fx_upper)), pch = 19, col = "blue", main = "Upper")
lines(x, pgamma(x, 2, 1))
## plot-lower
fx_lower <- discretize(pgamma(x, 2, 1),
method = "lower",
from = 0,
to = 10,
step = 0.5
)
x <- seq(0, 10, 0.5) #注意，lower 曲线绘制的时候只有三个参数
plot(stepfun(x, diffinv(fx_lower)), pch = 19, col = "blue", main = "Lower")
lines(x, pgamma(x, 2, 1))
### plot-rounding
fx_round <- discretize(pgamma(x, 2, 1),
method = "rounding",
from = 0,
to = 10,
step = 0.5
)
x <- seq(0, 10 - 0.5, 0.5)
plot(stepfun(x, diffinv(fx_round)), pch = 19, col = "blue", main = "Rounding")
lines(x, pgamma(x, 2, 1))
####plot-unbiased
fx_unbi <- discretize(pgamma(x, 2, 1),
method = "unbiased",
from = 0,
to = 10,
step = 0.5,
lev = levgamma(x, 2, 1) # Computes E[min(X, a)]
)
x <- seq(0, 10, 0.5)
plot(stepfun(x, diffinv(fx_unbi)), pch = 19, col = "blue", main = "Unbiased")
lines(x, pgamma(x, 2, 1))

###4.1.2 Calculation of the aggregate claim amount distribution

###总索赔金额分配的计算

#To perform Panjer’s recursion in R, we can use the function aggregateDist() with argument method
="recursive". 函数聚合+参数方法（递归）

#Example
We consider S such that N is Poisson distributed with mean 10 and X ∼ Gamma(2, 1). Then, we
approximate the cdf of S by first discretizing the gamma distribution on (0, 22) with the unbiased
method and a step of 0.5 and then using the recursive method in aggregate Dist()

fx <- discretize(pgamma(x, 2, 1),

method = "unbiased",
from = 0, to = 22, step = 0.5,
lev = levgamma(x, 2, 1)
)
Fs <- aggregateDist("recursive",
[Link] = "poisson",
[Link] = fx, lambda = 10,
[Link] = 0.5
)
Fs
#上面的代码返回一个类聚合目录的对象，我们可以从中获得关于 S 的信息

#use directly such an object to evaluate the cdf of S

Fs(20)

# the support can be obtained with the knots() function

head(knots(Fs))

#此外，我们还可以对上述对象应用其他函数，以获得不同的数量，甚至是图。
#summary() provides some basic information on S
summary(Fs)

#To plot the approximated cdf of S

plot(Fs, [Link] = FALSE, verticals = TRUE)

#Exact calculation by numerical convolutions 通过数值卷积进行的精确计算

#claim number 可以采用任何离散分布

#example
we consider a Gamma(2, 1) distribution for the severities, but now consider a Bin(10, 0.4)
distribution for the number of claims.

fx <- discretize(pgamma(x, 2, 1),

method = "unbiased",
from = 0, to = 22, step = 0.5,
lev = levgamma(x, 2, 1)
)
fn <- dbinom(0:10, 10, 0.4) # We require a vector of claim number probabilities
Fs <- aggregateDist("convolution",
[Link] = fn,
[Link] = fx,
[Link] = 0.5
)
Fs(10)

#density can be found using the diff() function

fs <- diff(Fs)
head(fs)

#Normal approximation
Fs <- aggregateDist("normal", moments = c(6, 2))
Fs(5)
# the above cdf evaluation can also be computed simply as
pnorm(5, mean = 6, sd = sqrt(2))

#Normal Power II approximation

Fs <- aggregateDist("npower", moments = c(6, 2, 0.5))
Fs(7) # Accesible only for x > mu_S

#Simulation
#The above can be performed using the aggregate Dist() function with the argument method =
"simulation".
#For instance, for S with N Poisson distributed with mean 10 and X ∼ Gamma(2, 1), we have
[Link](1)
model_freq <- expression(data = rpois(10))
model_sev <- expression(data = rgamma(2, 1))
Fs <- aggregateDist("simulation",
[Link] = 2500,
model_freq, model_sev
)

#We can then compute different quantities related to S. （平均值，分位数）

#For instance, the mean and quantiles can be computed using the mean() and quantile() functions:
mean(Fs)
quantile(Fs)

#Other relevant quantities in insurance applications are

#the Value at Risk (VaR) and the Conditional Tail Expectation.
VaR(Fs)

CTE(Fs)

#In fact, aggregateDist() implicitly calls the function simul() to perform the simulation.
#We can use the later directly as follows:
sim_s <- simul(2500,
[Link] = expression(rpois(10)),
[Link] = expression(rgamma(2, 1))
)

#This creates an object containing the severities, frequencies, and aggregate claim amounts.
(严重程度，频率，总索赔金额）
#The latter,representing S, can be obtained using the aggreate() function
s_sample <- aggregate(sim_s)
summary(s_sample[1, -1])
hist(s_sample[1, -1], freq = F, breaks = 20, main = "Simulation of S", xlab = "x")

#Alternatively, the simulation of S can be performed with the following code

s_sim <- function(n, distr, ...) {
sum(distr(n, ...))
}
[Link](1)
N <- rpois(2500, 10)
s_sample <- sapply(N, s_sim, distr = rgamma, shape = 2, rate = 1)
summary(s_sample)
hist(s_sample, freq = F, breaks = 20, main = "Simulation of S", xlab = "x")

###4.2 Ruin theory

###4.2.1 The surplus process
无 code
###4.2.2 The adjustment coefficient

在 R 中，我们可以使用 actuar package 中的 adjCoef（）函数来计算调整系数 adjustment coefficient

我们需要以下参数：两个时刻生成函数 MX (t)和 MW(t)（从而假设独立性），溢价 c，和上限的上限 MX (t)或任何其他上限 r。

#例如，如果 W∼Exp (2)，X∼Exp (1)和溢价是 c = 2.4，那么调整系数：

adjCoef(
 [Link] = mgfexp(x), [Link] = mgfexp(x, 2),
[Link] = 2.4, [Link] = 1
)

#In the above solution, we passed the upper bound for the support of MX(t).
#However, since we are in the classical risk model setup, we could have also used the upper bound
of R described previously
在上面的解中，我们通过了支持 MX (t)的上界。
然而，由于我们处于经典的风险模型设置中，我们也可以使用前面描述的 R 的上界

exp_aux <- function(x) {

x * dexp(x, 1)
}
snd_aux <- function(x) {
x^2 * dexp(x, 1)
}
exp_x <- integrate(exp_aux, 0, Inf)$value
exp_x

snd_x <- integrate(snd_aux, 0, Inf)$value

snd_x

c <- 2.4
lambda <- 2
bound <- 2 * (c - lambda * exp_x) / (lambda * snd_x)
bound

R <- adjCoef(
[Link] = mgfexp(x), [Link] = mgfexp(x, 2),
[Link] = 2.4, [Link] = bound
)
R

#knowledge of the adjustment coefficient allows computing a bound for the ruin probability,
# compute a bound for the ruin probability as (assuming u = 3)
u <- 3
exp(-R * u)

###4.2.3 Probability of ruin

#Let us illustrate first an exponential/exponential model with premium rate c = 1 (default):

psi <- ruin(
claims = "e", [Link] = list(rate = 5),
 wait = "e", [Link] = list(rate = 3)
)
psi(0:10) # Evaluates the ruin probability for initial surplus from 0 to 10.

#Next, we consider a model with a mixture of two exponentials for the claim amounts, exponential
interarrival times, and premium rate c = 1.5.
psi <- ruin(
claims = "e", [Link] = list(rate = c(3, 7), w = c(0.4, 0.6)),
wait = "e", [Link] = list(rate = 3),
pre = 1.5
)
psi(0:10)

Finally, we consider a model with Erlang claim amounts and exponentials interarrival times:
#we can plot the ruin probability straightforwardly as a function of the initial surplus
#using the plot() function as follows:
plot(psi, from = 0, to = 10)

#Finally, we consider a model with Erlang claim amounts and exponentials interarrival times:
psi <- ruin(
claims = "E", [Link] = list(shape = 3, rate = 1),
wait = "e", [Link] = list(rate = 3),
pre = 10
)
plot(psi, from = 0, to = 20)

###4.2.4 Reinsurance

#Proportional reinsurance 比例再保险

#Example 4.2.
Consider the following classical risk model under proportional reinsurance: The claim
amounts are exponentially distributed with mean 1, the Poisson rate is λ = 2, and the safety
loadings are θ = 0.2 and θh = 0.3.
a) Find the adjustment coefficient Rh if a = 0.75, 0.8, 0.9, 1.
b) Plot the adjustment coefficient as a function of the proportion a.
c) Find an upper bound for the ruin probability if a = 0.5 and u = 2.

#(a)
lambda <- 2
theta <- 0.2
thetah <- 0.3
# We require a function to compute the premium rate for different values of a
prem <- function(x) {
 ((1 + theta) - (1 + thetah) * (1 - x)) * lambda
}
# We require need a function to computute the mgf of aX for different values of a
mgf_ax <- function(x, y) {
mgfexp(x * y)
}
adj_prop <- adjCoef(mgf_ax,
[Link] = mgfexp(x, 2),
[Link] = prem,
upper = 1,
reins = "prop", from = 0, to = 1
)
adj_prop(c(0.75, 0.8, 0.9, 1))

#(b)
plot(adj_prop)

#(c)
exp(-adj_prop(0.5) * 2)

#Excess-of-loss reinsurance 超额损失再保险

Example 4.3.
Consider the following classical risk model under excess-of-loss reinsurance: The claim amounts
are Gamma distributed with shape parameter 2 and rate parameter 2, the Poisson rate is λ = 1,and
the safety loadings are θ = 0.2 and θh = 0.3. Plot the adjustment coefficient as a function of the
retention limit m varying from 0 to 10.

prem <- function(x) {

1.3 * levgamma(x, 2, 2) - 0.1
}
mgfx <- function(x, l) {
mgfgamma(x, 2, 2) * pgamma(l, 2, 2 - x) +
exp(x * l) * pgamma(l, 2, 2, lower = FALSE)
}
adj_eol <- adjCoef(mgfx,
premium = prem,
upper = 1,
reins = "excess-of-loss",
from = 0, to = 10
)
plot(adj_eol)