claim_model_3.R

claim_model_3.R

###########高速フーリエ変換によるクレーム総額の導出##########
library(actuar)

##FFTと一様分布畳み込み
#一様分布二つの和
x=discretize(punif(x),from=0, to=2,step=2^(-8),method="unbiased",lev=levunif(x))
y=fft(x)                           #密度関数ベクトルをフーリエ変換
z=y^2　　　　　　　　　　　　      #畳み込むため２乗する
fu=fft(z,inverse=TRUE)/length(z) 　#逆変換
fu=as.numeric(fu)　　　　　　　　　#虚数をとる
plot(seq(0, 2, by = 2^(-8)),fu,type="l",xlab="conv",ylab="frequency")

#一様分布三つの和
x=discretize(punif(x),from=0, to=3,step=2^(-8),method="unbiased",lev=levunif(x))
y=fft(x)                           #密度関数ベクトルをフーリエ変換
z=y^3　　　　　　　　　　　　      #畳み込むため３乗する
fu=fft(z,inverse=TRUE)/length(z) 　#逆変換
fu=as.numeric(fu)　　　　　　　　　#虚数をとる
plot(seq(0, 3, by = 2^(-8)),fu,type="l",xlab="conv",ylab="frequency")

##高速フーリエ変換でクレーム総額分布を求める
#クレーム額分布の確率密度の離散化
x=discretize(pgamma(x,1,0.1),from=0, to=512,step=0.0625,method="unbiased",lev=levgamma(x,1,0.1))
y=fft(x)                            #密度関数ベクトルをフーリエ変換
z=0.1/(1-0.9*y)                     #Nの確率密度関数に代入し、Sの特性関数ベクトルを得る
f8192=fft(z,inverse=TRUE)/length(z) #逆変換
f8192=as.numeric(f8192)             #虚数をとる
plot(seq(0, 512, by = 0.0625),f8192,type="l",xlab="S",ylab="frequency")
curve(0.009*exp(-0.01*x)*0.0625,0,512,add=TRUE,col=2)

##再帰式とFFTの比較
fx=discretize(pgamma(x,1,0.1),from=0, to=500,step=0.25,method="unbiased",lev=levgamma(x,1,0.1))
fx0=fx[1]

fs=rep(0,2001)

fs[1]=0.1
fs[2]=0.9*(fx[2]*fs[1])/(1-0.9*fx0)
fs[3]=0.9*(fx[3]*fs[1]+fx[2]*fs[2])/(1-0.9*fx0)

for(i in 2:2001){
k=rep(0,2001)
	for(j in 2:i){
	k[j]=fx[j]*fs[i-j+1]
	}
fs[i]=sum(k)*0.9/(1-0.9*fx0)
}


T=seq(0, 500, by = 0.25)
FFs=cumsum(fs)

plot(T,FFs,xlim=c(0,500),ylim=c(0,1),col=1,type="l")

par(new=TRUE)    

F8192=cumsum(f8192)
plot(seq(0, 500, by = 0.0625),F8192[1:8001],type="l",ylim=c(0:1),col=4)
curve(1-0.9*exp(-0.01*x),0,500,add=TRUE,col=2)