Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
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             #畳み込むため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             #畳み込むため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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment