ACTU-457 (Mast 724) Risk Theory

#example 
q <- 0.01
n1 <- 50
n2 <- 25
p.acc <- 0.3
B1.ord <- 50000
B1.acc <- 100000
B2.ord <- 75000
B2.acc <- 150000
#Group 1
mu.1 <- p.acc*B1.acc + (1-p.acc)*B1.ord
var.1 <- (p.acc*B1.acc^2+(1-p.acc)*B1.ord^2)-mu.1^2
#Group 2
mu.2 <- p.acc*B2.acc + (1-p.acc)*B2.ord
var.2 <- (p.acc*B2.acc^2+(1-p.acc)*B2.ord^2)-mu.2^2



#
x <- c(0,1,2,3,4,5,6,7,8,9,10)
px <- c(0,0.15,0.2,0.25,0.125,0.075,0.05,0.05,0.05,0.025,0.025)
n <- c(0,1,2,3,4,5,6,7,8)
pn <- c(0.05,0.1,0.15,0.2,0.25,0.15,0.06,0.03,0.01)
EX <- sum(x*px)
EN <- sum(n*pn)
ES <- EX*EN*25
EX2 <- sum(x^2*px)
EN2 <- sum(n^2*pn)
VarX <- EX2-EX^2
VarN <- EN2-EN^2
25*sqrt(EN*VarX+VarN*EX^2)

[1] 191.6155

cov.ite.prec <- px
size.convo <- 11+10*n
size.max <- 11+10*9
conv.ite <- matrix(c(1,rep(0,21),rep(0,size.max-22),px,rep(0,11),rep(0,size.max*10-44-79)),nrow=size.max,ncol=10)
convolution <- px 


for (i in 1:8){
    conv.ite[,i+2] <- round(c(convolve(convolution,rev(px),type="open"), rep(0,size.max-length(convolve(convolution,rev(px),type="open")))),5)
    convolution <- conv.ite[1:(10+i*9),i+2] 
}
convo.table <- conv.ite[1:22,1:9]
convo.table

      [,1]  [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]    [,9]
 [1,]    1 0.000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
 [2,]    0 0.150 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
 [3,]    0 0.200 0.02250 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
 [4,]    0 0.250 0.06000 0.00338 0.00000 0.00000 0.00000 0.00000 0.00000
 [5,]    0 0.125 0.11500 0.01350 0.00051 0.00000 0.00000 0.00000 0.00000
 [6,]    0 0.075 0.13750 0.03488 0.00270 0.00008 0.00000 0.00000 0.00000
 [7,]    0 0.050 0.13500 0.06144 0.00878 0.00051 0.00001 0.00000 0.00000
 [8,]    0 0.050 0.10750 0.08569 0.01999 0.00198 0.00009 0.00000 0.00000
 [9,]    0 0.050 0.08812 0.09750 0.03580 0.00549 0.00042 0.00002 0.00000
[10,]    0 0.025 0.07875 0.09841 0.05266 0.01194 0.00136 0.00008 0.00000
[11,]    0 0.025 0.07062 0.09337 0.06682 0.02138 0.00345 0.00031 0.00002
[12,]    0 0.000 0.06250 0.08812 0.07597 0.03282 0.00726 0.00091 0.00007
[13,]    0 0.000 0.04500 0.08370 0.08068 0.04450 0.01305 0.00218 0.00022
[14,]    0 0.000 0.03125 0.07673 0.08266 0.05485 0.02062 0.00448 0.00060
[15,]    0 0.000 0.01750 0.06689 0.08277 0.06314 0.02930 0.00808 0.00138
[16,]    0 0.000 0.01125 0.05376 0.08081 0.06933 0.03826 0.01304 0.00280
[17,]    0 0.000 0.00750 0.04125 0.07584 0.07361 0.04677 0.01919 0.00505
[18,]    0 0.000 0.00500 0.03052 0.06811 0.07578 0.05438 0.02616 0.00830
[19,]    0 0.000 0.00313 0.02267 0.05854 0.07552 0.06080 0.03352 0.01254
[20,]    0 0.000 0.00125 0.01674 0.04878 0.07262 0.06573 0.04084 0.01768
[21,]    0 0.000 0.00063 0.01167 0.03977 0.06747 0.06882 0.04775 0.02351
[22,]    0 0.000 0.00000 0.00766 0.03184 0.06079 0.06982 0.05389 0.02978

fs <- 
function(s){
sum(pn[1:9]*convo.table[s+1,1:9])
}
fs(1)

[1] 0.015

fs.vec<- pn[1]
for(i in 1:21){fs.vec[i+1] <- fs(i)}

FS.vec <- pn[1]
FS.vec <- cumsum(fs.vec[1:22])

ES <- EX*EN*25
LimExp25 <- ES-25*(1-FS.vec[2])
LimExp50 <- LimExp25-25*(1-FS.vec[3])



#fs525 <- pn[1]*1 + 
#       pn[2]*px[1]+
#       pn[3]*conv.ite[1,1]+pn[2]*px[2] +
#       ...
#fs525 = 0.2479
#one may use sumprod    


314.5-25*(1-0.065)

[1] 291.125

291.125-25*(1-0.08838)

[1] 268.3345

#Example 
p1 <- (1-0.4)*dbinom(1,3,0.3)/(1-dbinom(0,3,0.3))
p0 <- dbinom(0,3,0.3)
a <- -3/7
b <- 12/7
fX <- c(0.3,0.5,0,0.2)
fS <-   0.53702
fS1 <- ((p1-(a+b)*0.4)*fX[2]+(a+b*1/1)*fX[2]*fS)/(1-a*0.3)
fS2 <- ((p1-(a+b)*0.4)*fX[3]+
                (a+b*1/2)*fX[2]*fS1+
                (a+b*2/2)*fX[3]*fS)/(1-a*0.3)

EY <- (1/gamma(2))*sqrt(10000)*gamma(1.5)*gamma(1.5)
ER <- 1*(0.5)+1.3*(0.4)+2*0.1
EX <- EY*ER

BurrDist <- function(y,lam=10000,aa=2,tt=2){return(1-(lam/(lam+y^tt))^aa)}
0.5*BurrDist(200)+0.4*BurrDist(153.8462)+0.1*BurrDist(100)

[1] 0.9197135

###################################################################
EBurr <- function(x,aa=2,tt=2,lam=10000){return(x*aa*tt*lam^aa*x^(tt-1)/(lam+x^tt)^(aa+1))}
#verification EY
integrate(EBurr,0,100000)$value

[1] 78.53982

###########################
ELbound200 <- integrate(EBurr,200,100000)$value
ELbound153 <- integrate(EBurr,153.8462,100000)$value
ELbound100 <- integrate(EBurr,100,100000)$value
0.5*(ELbound200-200*(1-BurrDist(200)))+
0.4*1.3*(ELbound153-153.8462*(1-BurrDist(153.8462)))+
0.1*2*(ELbound100-100*(1-BurrDist(100)))

[1] 7.550433