1. 假设某险种的损失记录如下。
年底 损失次数 平均损失金额（元）
2017 200 1200
2018 400 1500

——2380.2

  1. 假设某保单的损失服从指数分布，密度函数为 \[f(x;\lambda)=\lambda\exp(-\lambda x), x>0,\] 其中 \(\lambda\) 为未知参数。如果该保单过去各年的损失观测值为 \(x_1, x_2, ..., x_n\)，求参数 \(\lambda\) 的极大似然估计值。

——-1/E(x)

  1. 计算混合指数分布 \[f(x)=\sum_{i=1}^{n}a_i \lambda_i \exp(-\lambda_i x), x>0, \lambda_i >0, a_i >0, \sum_{i=1}^{n}a_i=1\] 的矩母函数。
  1. 假设保单规定的免赔额为 20，而保单的损失服从参数为 0.2 的指数分布，求保险人对该保单的期望赔款。
  1. 假设一组损失数据的均值为 100，标准差为 223.607。
#正态分布
Var.nor<-function(q)qnorm(q,mean=100,sd=223.607)
TVar.nor<-function(q){
  integrate(Var.nor,lower=q,upper=1)$value/(1-q)
}
Var.nor(0.9);TVar.nor(0.9)
## [1] 386.5639
## [1] 492.4266
Var.nor(0.99);TVar.nor(0.99)
## [1] 620.1877
## [1] 695.9606
Var.nor(0.999);TVar.nor(0.999)
## [1] 790.9976
## [1] 852.9049
#帕累托分布
library(actuar)
## Warning: package 'actuar' was built under R version 3.6.3
## 
## Attaching package: 'actuar'
## The following object is masked from 'package:grDevices':
## 
##     cm
Var.par<-function(q) qpareto(q,shape=2.5, scale=150)
TVar.par<-function(q){
  integrate(Var.par,lower=q,upper=1)$value/(1-q)
}
Var.par(0.9);TVar.par(0.9)
## [1] 226.783
## [1] 477.9716
Var.par(0.99);TVar.par(0.99)
## [1] 796.436
## [1] 1427.393
Var.par(0.999);TVar.par(0.999)
## [1] 2227.34
## [1] 3812.233
#威布尔分布
qwei<-function(y,alpha,theta){
  f<-((-log(1-y))/alpha)^(1/theta)
  return(f)
}
Var.wei<-function(q) qwei(q,alpha=0.1*sqrt(2),theta=0.5)
TVar.wei<-function(q){
  integrate(Var.wei,lower=q,upper=1)$value/(1-q)
}
Var.wei(0.9);TVar.wei(0.9)
## [1] 265.0949
## [1] 595.3534
Var.wei(0.99);TVar.wei(0.99)
## [1] 1060.38
## [1] 1620.897
Var.wei(0.999);TVar.wei(0.999)
## [1] 2385.854
## [1] 3176.63
6. 假设被保险人的损失 \(X\) 服从伽马分布，其中形状参数为 2， 尺度参数为 1000。两份保单如下： - 保单 A 的免赔额为100。 - 保单 B 的免赔额为100，赔偿限额为 3000。（\(d=100, u=3100\)) - 分别计算保险公司对保单 A 和保单 B 的期望赔款（含零赔款在内）。 - 如果发生10% 的通货膨胀，上述结果将如何变化？ - 如果通胀函数为 \(1.1x^{0.5}\)，上述结果将如何变化？ 
#伽马分布的生存函数
S<-function(x) 1-pgamma(x,2,1/1000)
#伽马分布的期望赔款（含零赔款在内）
#（1）
#A--只含免赔额100
A1<-integrate(S,100,Inf)$value
A1
## [1] 1900.159
#B--免赔额100，限额3000
B1<-integrate(S,100,3100)$value
B1
## [1] 1670.408
# (2)
# A
A2<-(integrate(S,100/1.1,Inf)$value)*1.1
A2
## [1] 2100.132
# B
B2<-(integrate(S,100/1.1,3100/1.1)$value)*1.1
B2
## [1] 1783.645
# (3)
S<-function(x) 1-pgamma(x,2,1/1000)
#A--只含免赔额100
A3<-(integrate(S,(100/1.1)^2,Inf)$value)^0.5*1.1
A3
## [1] 1.788362
#B--免赔额100，限额3000
B3<-(integrate(S,(100/1.1)^2,(3100/1.1)^2)$value)^0.5*1.1
B3
## [1] 1.788362
7. 对于伽马分布，其中形状参数为 2， 尺度参数为 1000，绘制以下图形： - 止损保费和平均超额损失随着免赔额增加而变化的曲线图 - 有限期望值随着限额变化而变化的曲线图 - 把上述分布改为 帕累托\((\alpha = 2, \theta = 200)\) 和 指数分布 (\(\theta = 1/200\)) - 注意： 可以使用 actuar 包的 ppareto2 函数。 
#（1） 伽马分布（形状参数为 2， 尺度参数为 1000）
library(actuar)
#伽马分布的生存函数
S<-function(x) 1-pgamma(x,2,1/1000)
#伽马分布的止损保费函数 x1
x1<-NULL
d1<-seq(1,150,0.1)
for(i in 1:length(d1)){
  x1[i]<-integrate(S,d1[i],Inf)$value
}
#伽马分布的平均超额函数 ex1
ex1<-NULL
for(i in 1:length(d1)){
  ex1[i]<-integrate(S,d1[i],Inf)$value/S(d1[i])
}
#伽马分布的有限期望值 m1
S<-function(x) 1-pgamma(x,2,1/1000)
m1<-NULL
for(i in 1:length(d1)){
  m1[i]<-integrate(S,0,d1[i])$value
}
#绘图
par(mfrow=c(1,2))
plot(d1,x1,type='l',main='伽马分布',xlab='x',ylab='y')
lines(d1,ex1,type='l',col=2,lty=2)
legend('topright',c('止损保费','平均超额损失'),lty=c(1,2),col=c(1,2))
plot(d1,m1,type='l',main='伽马分布',xlab='x',ylab='y')
legend('topleft','有限期望值',lty=1,col=1)

#指数分布（theta = 1/200）
S<-function(x) exp(-x*(1/200))
#指数分布的止损保费函数 x1
x1<-NULL
d1<-seq(0.1,150,0.1)
for(i in 1:length(d1)){
  x1[i]<-integrate(S,d1[i],Inf)$value
}
#指数分布的平均超额函数 ex1
ex1<-NULL
for(i in 1:length(d1)){
  ex1[i]<-integrate(S,d1[i],Inf)$value/S(d1[i])
}

#指数分布的有限期望值 m1
S<-function(x) exp(-x*(1/200))
m1<-NULL
for(i in 1:length(d1)){
  m1[i]<-integrate(S,0,d1[i])$value
}
#绘图
plot(d1,x1,type='l',main='指数分布',xlab='x',ylab='y')
lines(d1,ex1,type='l',col=2,lty=2)
legend('bottomright',c('止损保费','平均超额损失'),lty=c(1,2),col=c(1,2))
plot(d1,m1,type='l',main='指数分布',xlab='x',ylab='y')
legend('topleft','有限期望值',lty=1,col=1)

#(3) 帕累托分布（alpha = 2, theta = 200）
#帕累托函数分布的生存函数
library(actuar)
alpha<-2
theta<-200  
S<-function(x){
 1-ppareto2(x,shape=2,scale=200)
}
#帕累托分布的止损保费函数 x1
x1<-NULL
d1<-seq(0.1,500,1)
for(i in 1:length(d1)){
  x1[i]<-integrate(S,d1[i],Inf)$value
}
#帕累托分布的平均超额函数 ex1
ex1<-NULL
for(i in 1:length(d1)){
  ex1[i]<-integrate(S,d1[i],Inf)$value/S(d1[i])
}
#帕累托分布的有限期望值 m1
m1<-NULL
for(i in 1:length(d1)){
  m1[i]<-integrate(S,0,d1[i])$value
}
#绘图
plot(d1,x1,type='l',main = '帕累托分布',xlab='x',ylab='y')
lines(d1,ex1,type='l',col=2,lty=2)
legend('topleft',c('止损保费','平均超额损失'),lty=c(1,2),col=c(1,2))
plot(d1,m1,type='l',main='帕累托分布',xlab='x',ylab='y')
legend('topleft','有限期望值',lty=1,col=1)

  1. 假设个体风险服从参数为 \(\left( \lambda \Theta \right)\) 泊松分布， 结构函数 \(\Theta\) 为均值为1，方差为 \(1/\alpha\) 伽马分布，相关参数表示为 \(\left( \alpha ,\alpha \right)\)，且密度函数为

\[{{f}_{\Theta }}(\theta )=\frac{1}{\Gamma (\alpha )}{{\alpha }^{\alpha }}{{\theta }^{\alpha -1}}\exp (-\alpha \theta ),\theta >0\]

alpha =1
lambda =1.2
x=seq(1,10,1)
dnbinom(x,size=1,mu=1.2)  
##  [1] 0.247933884 0.135236664 0.073765453 0.040235702 0.021946746 0.011970953
##  [7] 0.006529610 0.003561606 0.001942694 0.001059651
  1. 假设某保险业务的累积损失 \(S\) 服从复合泊松分布，泊松分布的参数为 20。每次损失金额服从均值为 100 的指数分布。用平移伽马近似求累积损失的 99% 的分位数。
x0 <- -666.67
alpha <- 17.778
beta <- 6.667*10^(-3)
x<-qgamma(0.99,shape=alpha, scale=1/beta)+x0
x
## [1] 3687.061
  1. 某保险公司承保了 160 个建筑物火灾保险，有关数据如下表。假设对一个建筑物，其发生火灾的概率是 0.04，每个建筑物的火灾为相互独立事件，在发生赔款的条件下，保险公司的赔款额是从 0 到最高赔款额的均匀分布。
建筑物种类 最高赔款额（元） 保单数
1 10000 80
2 20000 35
3 30000 25
4 50000 15
5 100000 5

#均值为710^4,方差为4.22410^4

  1. 损失次数 \(N\): 泊松分布 (\(\lambda = 3\))
# 对pareto离散化后，求复合分布S
dpareto <- function(x) {
  L = 4 * (10^4)/((x + 10)^(4 + 1))
  return(L)
}
ppareto <- function(x) {
  L = 1 - (10/(x + 10))^4
  return(L)
}
v = 1 - ppareto(6)
lambda = 3 * v
Yp <- function(y) {
  if (y < 13.5) 
    L = 1 - (1 - ppareto(y/0.75 + 6))/v else L = 1
    return(L)}
# pareto分布的离散化
h = 2.25
f = Yp(h/2)
k = 0
repeat {
  k = k + 1
  last = Yp((k + 0.5) * h) - Yp((k - 0.5) * h)
  f = c(f, last)
  if (sum(f) >= 0.99) 
    break
}
# 递推公式求S
g = exp(-lambda * (1 - f[1]))
k = 0
repeat {
  k = k + 1
  summ = 0
  if (k < length(f)) 
    upper = k else upper = length(f) - 1
  for (j in 1:upper) {
    summ = summ + j * f[j + 1] * g[k - j + 1]
  }
  g = c(g, summ * lambda/k)
  if (sum(g) > 0.999999) 
    break
}
sum(g)
## [1] 0.9999995
plot(g,type='h')

12. 假设损失次数服从负二项分布, 参数为 (\(r = 2, \beta = 3\)), - 每次损失的金额服从对数正态分布, 参数为 (\(\mu=5, \sigma=2\)), - 计算累积损失在 90%、95% 和 99% 水平下的 VaR 和 TVaR。 - 注：累积损失的分布用随机模拟。 
set.seed(321)
r = 2;beta = 3
mu = 5; sigma = 2     #对数正态分布的参数
s = n = NULL
for ( i in 1:10000) {
  n[i] = rnbinom(1,size=r, mu=r*beta)
  s[i] = sum(rlnorm(n[i], meanlog = mu, sdlog = sigma))
}
#绘图
s = sort(s)
hist(s, freq = F, breaks = 100, col = 2)

plot(s, cumsum(s)/sum(s), type = 's', col = 2)

VaR0.9=quantile(s,0.9)
TVR0.9=mean(s[s>VaR0.9])
VaR0.9
##      90% 
## 15097.24
TVR0.9
## [1] 40931.09
VaR0.95=quantile(s,0.95)
TVR0.95=mean(s[s>VaR0.95])
VaR0.95
##      95% 
## 25545.26
TVR0.95
## [1] 62518.93
VaR0.99=quantile(s,0.99)
TVR0.99=mean(s[s>VaR0.99])
VaR0.99
##     99% 
## 69786.8
TVR0.99
## [1] 158768.1
  1. 某团体意外伤害保险的有关数据如下：
类别 人数 概率 保险金（万元）
1 100 0.01 50
2 120 0.02 40
3 200 0.025 45
4 400 0.018 55
5 260 0.015 60
6 150 0.01 65
  1. 损失次数服从泊松 (\(\lambda = 3\))，
#随机模拟法
set.seed(111)
lam = 3                   #泊松参数
alpha = 4; theta =10       #帕累托分布的参数
u = 24
d=6
s = n = NULL                    # 保险人的年度累积赔款
# 开始模拟
for (i in 1:10000){
  n[i] = rpois(1, lambda = lam)  # 模拟损失次数
  x <- rpareto2(n[i], shape = alpha, scale = theta)
  w <- 0.75*pmax(x-d, 0)  # 保险人对每次损失的所负金额
  s[i] = sum(pmin(w, 0.75*(u-d)))
  s[i] = round(s[i])
}
hist(s, freq = F, breaks = 1000, col = 2)

s = sort(s)
plot(s, cumsum(s)/sum(s), type = 's', col = 2)