用r-也能做精算actuar包學(xué)習(xí)筆記

引 數(shù)據(jù)描 構(gòu)造分組數(shù)據(jù)對(duì) 分組數(shù)據(jù)分布 計(jì)算數(shù)據(jù)經(jīng)驗(yàn) 損失分 損失分布種 損失分布的估 損失隨量的修 風(fēng)險(xiǎn)理 復(fù)合分 連續(xù)分布的離散 復(fù)合分布的計(jì) VaR和 保單組合的模 復(fù)合層次模型的模 模擬結(jié)果的處 信度理 信度理論簡(jiǎn) 層次信度模 信度回歸模 ChapterChapter引對(duì)R學(xué)包actar單actuar于20052006年2開(kāi)入R1.1actuarR包中的很函數(shù)可以算師使用，是為了達(dá)到個(gè)目的而尋某個(gè)包的某函數(shù)是一個(gè)費(fèi)actuar如題所示，本文是我在學(xué)習(xí)ar程中的學(xué)習(xí)筆記，主要涉及這個(gè)包中一些函數(shù)的使而可能講的會(huì)比較詳細(xì)。文章主要是針對(duì)R語(yǔ)言的初學(xué)者，因此每種函數(shù)或數(shù)據(jù)的結(jié)構(gòu)進(jìn)行了盡可能直，正式開(kāi)數(shù)據(jù)構(gòu)造分組數(shù)據(jù)對(duì)損失數(shù)據(jù)的型主要分為組數(shù)據(jù)和組數(shù)據(jù)。對(duì)于組數(shù)據(jù)描述方法大會(huì)比有的處理原始數(shù)據(jù)的基本方法，通過(guò)將數(shù)據(jù)分到不同的組中，我們可以看到各組中數(shù)據(jù)的相對(duì)頻數(shù)，有助于對(duì)數(shù)據(jù)形成直觀的印象比如我們對(duì)連續(xù)變量繪制直方圖；而且在生存函數(shù)的估計(jì)中，再使用分組數(shù)據(jù)進(jìn)行生存函數(shù)的估計(jì)，這樣可以有效減小計(jì)算量?，F(xiàn)在假設(shè)我們要把一組連續(xù)變量分為rc0,c1,c1,c2,...,cr?1,cr，那么就需要定義r+1個(gè)邊界c0c1crc00造一個(gè)分組數(shù)據(jù)的結(jié)構(gòu)，而非對(duì)現(xiàn)有連續(xù)數(shù)據(jù)進(jìn)行分組，該函數(shù)返回一個(gè)分組數(shù)據(jù)的對(duì)象(groupeddatabject)。函數(shù)語(yǔ)法使用說(shuō)明Group定義的是分組的邊界值，freq1和freq2(還可以定義freq3及)是每條分組數(shù)據(jù)的頻數(shù)。Group，freq1freq2可以隨意命名，比如我們可以將Group改為“分?jǐn)?shù)檔”，將freq1，freq2Group向量要比f(wàn)req向量多出一個(gè)長(zhǎng)度(邊界數(shù)比組數(shù)多1)返回的是一個(gè)數(shù)據(jù)框。特別要注意對(duì)第一列的處理，見(jiàn)下面例例子options(digits=x=grouped.data(Group=c(0,25,50,100,150,250,500),Line.1= 31,57,42,65,84),Line.2=c(26,33,31,19,16,123456改成左閉右開(kāi)區(qū)間，并自定義行名稱(chēng)x1=grouped.data(Group=c(0,25,50,100,150,250,500),Line.1= 31,57,42,65,84),Line.2=c(26,33,31,19,16, right=F,s= BCDEF為避免修改原始數(shù)據(jù)x，以下x賦值到x2，對(duì)x2進(jìn)行操作x2=x2[1,GroupLine.1Line.21 提取第一條分組數(shù)據(jù)(某列)x2[,[1]3031574265提取各組的邊界值。如果第一列你期待會(huì)出現(xiàn)什么x2[, 50100150250的修改，一定要同時(shí)指定分組區(qū)間的左右的邊界值。比如下面這條命令將第一組的右邊界由2：(x2[1,1]=c(0, 0體會(huì)這樣修改波及的范圍(x2[c(3,4),1]=c(55,110, 551101。(x2[1,1]=[1]數(shù)據(jù)落在每組中的頻數(shù)呢？答案就是使用cut函數(shù)。例子1005z=rexp(100,rate=指定邊界點(diǎn)，也是劃分點(diǎn)break.points=c(0,1,4,8,14,(tz=table(cut(z,breaks=(8,14] 用匯總結(jié)果直接構(gòu)造分組數(shù)據(jù)對(duì)象grouped.data(Group=break.points,freq=Group123455分組數(shù)據(jù)分布有了grouped.data對(duì)象，我們就可以對(duì)該對(duì)象進(jìn)行一系列操作。首先是繪制分組數(shù)據(jù)的經(jīng)R會(huì)應(yīng)用分組數(shù)據(jù)對(duì)象x的第一列劃分繪圖時(shí)需要指定頻率所在的列，如果不指定，默認(rèn)繪制第一組頻率。例子layout(matrix(1:3,1,hist(x[,-3],main="Histogramofhist(x[,-2],main="Histogramofhist(x,main="HistogramforUnspecified0.0000.0010.0020.003 0.0000.0010.0020.0030.0000.0010.0020.003 0.0000.0010.0020.003 x[, x[, 繪制拱形圖。如同對(duì)連續(xù)的隨量可以繪制經(jīng)驗(yàn)分布函數(shù)圖一樣，對(duì)于分組數(shù)據(jù)也可(ogi x≤ (cj?x)Fn(cj?1)+(x?cj?1)cj

cj?1<x≤ x>函數(shù)e()輸入的是分組數(shù)據(jù)對(duì)象，返回的是一個(gè)階梯函數(shù)對(duì)象(Stepunctnlass)，也就是說(shuō)實(shí)現(xiàn)了分組數(shù)據(jù)對(duì)象向階梯函數(shù)對(duì)象的轉(zhuǎn)換。如果給定函數(shù)的橫坐標(biāo)，就可以ecdfots/t例子得到一個(gè)階梯函數(shù)Fnt=返回臨界點(diǎn) 50100150250返回臨界點(diǎn)對(duì)應(yīng)的累積頻率值[1]0.000000.097090.197410.381880.517800.72816對(duì)函數(shù)作圖，得到ogive●●●●●●●●●●●●● x計(jì)算數(shù)據(jù)經(jīng)驗(yàn)首先是計(jì)算經(jīng)驗(yàn)1一階矩。函數(shù)mean是一個(gè)泛型函數(shù)(genericfunction)2，除可以作用于通常的組數(shù)據(jù)向量對(duì)象(vector)以外，還可以作用于分組數(shù)據(jù)對(duì)象(grouped.data)，計(jì)算 1n

cj?1+2

r個(gè)組的總值相加再除以樣本量n，就得到每個(gè)觀測(cè)均值的估計(jì)。該公式假設(shè)每組內(nèi)的觀測(cè)值分布是均勻的。例子以上文中的分組數(shù)據(jù)對(duì)象x為例x1231文中經(jīng)常會(huì)提到“經(jīng)驗(yàn)**”，比如經(jīng)驗(yàn)分布函數(shù)，經(jīng)驗(yàn)一階矩等，個(gè)人理解所謂經(jīng)驗(yàn)就是將樣本當(dāng)成總體去對(duì)待，比如經(jīng)驗(yàn)方差是除以樣n，而樣本方差是n-1bootstrap方法中，使用經(jīng)驗(yàn)分布函數(shù)去代替總體2又譯作類(lèi)函數(shù)，泛型函數(shù)將作用對(duì)象的屬性也作為一個(gè)參數(shù)輸入，對(duì)于不同的對(duì)象類(lèi)別采用不同的方法，并得到同的輸出。actuar包中使得mean函數(shù)可以作aggregateDistgrouped.data對(duì)象，在加載actuar包后，可以通過(guò)命令methods(mean)查看mean函數(shù)的所有方法。456Line.1 ((0+25)/2*30+(25+50)/2*31+(50+100)/2*57+(100 150)/2*42++(150+250)/2*65+(250+500)/2* [1]mean()emm函數(shù)則可以計(jì)算任意階的經(jīng)驗(yàn)原點(diǎn)矩。首先引入actuar包中的兩個(gè)數(shù)據(jù)集。其中dental是組數(shù)據(jù)，gdental是分組數(shù)據(jù)。 403512593171511107cj 25] ( 50] (50,100] (100,150] (150,250] (250,500] 89103 其中，order是階數(shù)，可以賦值給它一個(gè)向量，這樣就能計(jì)算多個(gè)原點(diǎn)矩。x可以是數(shù)據(jù)向量或者是矩陣，對(duì)于矩陣，emm將每一列視為一條數(shù)據(jù)。組數(shù)據(jù)k階經(jīng)驗(yàn)矩的計(jì)算公式為例子

組數(shù)據(jù)向量形式

1n

xjxemm(dental,[1]數(shù)據(jù)矩陣形xx=matrix(1:9,3,147258369xxxx第二列的均值和emm(xx,258如果是分組數(shù)據(jù)，xgrouped.data()k階經(jīng)驗(yàn) 例子emm(gdental,

1∑n

?cj?1j(k+1)·(cj?

[1]3.533e+023.577e+05u(2.3節(jié))u的損失額度被強(qiáng)制定義為u。elev函數(shù)可以計(jì)算經(jīng)驗(yàn)有限期望值(empiricallimitedexpectedvalue)，其中x可以是組數(shù)據(jù)，也可以是分組數(shù)據(jù)。1n

min(xj, 在此略去，有的同學(xué)可以參考幫助文檔。uu計(jì)算出的經(jīng)驗(yàn)有限期望值是不一樣的。elevu的經(jīng)驗(yàn)有限期望值，只需要特別指定u即可。例子組數(shù)據(jù)，返回的lev是上限u的函數(shù)lev=這里將保單限額u200[1]lev函數(shù)的拐點(diǎn)，注意到拐點(diǎn)都發(fā)生在數(shù)據(jù)點(diǎn)，而觀察下圖可知拐點(diǎn)間的函數(shù)都是線 46107141259317351567分組數(shù)lev2=[1]分別對(duì)組數(shù)據(jù)和分組數(shù)據(jù)的有限期望函數(shù)作圖par(mfrow=c(1,plot(lev,type="o",pch=plot(lev2,type="o",pch=●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●Empirical200Empirical100 100 Chapter損失損失分布種根據(jù)損失額的特征，損失分布常選用具有非負(fù)支集(密度函數(shù)f(x)的支集指的是使得f(x)?=0x的集合)的連續(xù)分布。Rd，p，q，r四種函數(shù)，分別是密度函數(shù)、分布函數(shù)、分布函數(shù)的反函數(shù)(分位數(shù))和生成該分布的隨機(jī)數(shù)。actuar包提供了與[5]的附錄A中所列示的連續(xù)分布族相配套的這四種函數(shù)(除去逆和對(duì)數(shù)t分布，但包括Gamma分布)Rstats中并不自帶，但有些分布在精算研究中卻很(pareto分布)。此外，actuarm、lev和mgf三種函數(shù)，m是計(jì)算理論原點(diǎn)矩，lev是計(jì)算有限期望值，mgf是計(jì)算矩母函數(shù)。密度函數(shù)、分布函數(shù)、原點(diǎn)矩、有限期望值及其k次方都可以通過(guò)查詢(xún)?cè)摳戒浀玫健?duì)于經(jīng)驗(yàn)數(shù)據(jù)，如上面所介紹的，actuaremmelev來(lái)計(jì)算經(jīng)驗(yàn)原點(diǎn)矩和經(jīng)驗(yàn)有限期望值(這兩個(gè)函數(shù)的前綴都是empirical)。ratescale參數(shù)，scale=1/rate，因此兩者在本質(zhì)上是等價(jià)的，[5]中使用的是scale參數(shù)，在指定參數(shù)時(shí)千萬(wàn)不要弄混。例子這里以雙參數(shù)pareto分布為par(mfrow=c(1,繪制密度函數(shù)曲線curve(dpareto(x,shape=2,scale=2),from=0.001,to= main="density繪制分布函數(shù)曲線curve(ppareto(x,shape=2,scale=2),from=0.001,to= main="cumulativedistributiondpareto(x,shape=2,scale=ppareto(x,dpareto(x,shape=2,scale=ppareto(x,shape=2,scale= pareto分布中位數(shù)qpareto(0.5,shape=2,scale=[1]5paretorpareto(5,shape=2,scale=[1]7.910250.098170.188240.48281E(X1.5)，注pareto分布是厚尾分布k階矩要求?1<k<ααshape參mpareto(order=1.5,shape=2,scale=[1]E[(X5)1.5]，注意同?1kαlevpareto(limit=5,shape=2,scale=2,order=[1]curve(mgfexp(x,rate=2),-1,mgfexp(x,rate=0.81.0mgfexp(x,rate=0.81.01.8x損失分布的估R中，MASS?tdistr函數(shù)可以進(jìn)行極大似然估計(jì)。在actuar包中，mde函數(shù)則提供了三種距離最小化的分布擬合方法(miniumdistanceestimates)。Cram′er-vonMises(CvM)(對(duì)于分組數(shù)據(jù)是ogive)的距離。未分組數(shù)據(jù)d(θ)

wj(F(xj;θ)? 分組數(shù)據(jù)

d(θ)

c 在這里，F(xiàn)(x;θ)分布函數(shù)，θ是其參數(shù)；Fn(x)是經(jīng)驗(yàn)分布函數(shù)ecdf；F?n(x)是分組數(shù)據(jù)的經(jīng)驗(yàn)分布函數(shù)ogive；wj是賦予每個(gè)觀測(cè)或組別的權(quán)重，默認(rèn)都取1。修正卡方法僅應(yīng)用于分組數(shù)據(jù)，通過(guò)最小化各組期望頻數(shù)與實(shí)際觀測(cè)頻數(shù)的平方誤差到d(θ)

wj[n(F(cj;θ)?F(cj?1;θ))?j 其中n

njwj默認(rèn)情n?1jLAS(layeraverageseverity)也僅應(yīng)用于分組數(shù)據(jù)。通過(guò)最小化各組內(nèi)的理論和經(jīng)驗(yàn)jd(θ)

wj(LAS(cj?1,cj;θ)?LA?Sn(cj?1, L?Sn(x,?n?njjstatsoptim函數(shù)語(yǔ)法mde(x,fun,start,measure=c("CvM","chi-使用說(shuō)明x是分組數(shù)據(jù)對(duì)象的或未分組的數(shù)據(jù)fun是待擬合的分布，CvM法和修正卡方法需要指定分布函數(shù)：p**。LAS法需要指定理論有限期望函數(shù)lev**。start指定參數(shù)初始值。形式必須以列表的形式，形式可以見(jiàn)例子，有幾個(gè)參數(shù)就要指定measure是指定方法。weight指定權(quán)重，否則采用默認(rèn)權(quán)optimL-BFGS-B方法進(jìn)行優(yōu)化可以添加參數(shù)method=“L-BFGS-B”。mde(list)，rate是參數(shù)估計(jì)結(jié)果，distance是最小化后的距離。我們可以對(duì)上面的gdental數(shù)據(jù)進(jìn)行分布擬合，這是一個(gè)分組數(shù)據(jù)，在為參數(shù)估計(jì)賦初始值時(shí)，假設(shè)數(shù)據(jù)來(lái)自均值為200的指數(shù)分布。。例子首先觀察一下數(shù)據(jù)的分布CvM(mde.est1=mde(gdental,pexp,start=list(rate=1/200),measure=mu1=修正卡方法(mde.est2=mde(gdental,pexp,start=list(rate=1/200),measure="chi-mu2=LAS(mde.est3=mde(gdental,levexp,start=list(rate= measure=mu3=CvMLAScurve(mu1*exp(1)^(-mu1*x),from=0,to=4000,add=+col=curve(mu2*exp(1)^(-mu2*x),from=0,to=4000,add=+col=curve(mu3*exp(1)^(-mu3*x),from=0,to=4000,add=+col=legend(2700,0.0025,legend=c("CvM","chi-square",+col=c("red","blue","green"),lty= 我們還可以對(duì)組數(shù)據(jù)進(jìn)行分布擬合，下面的一個(gè)例子是一個(gè)混合分布例子首先生成400個(gè)隨機(jī)數(shù)，其中200個(gè)來(lái)自Gamma(α2θ2)，200個(gè)來(lái)自Gamma(α=10,θ=2)。dat=c(rgamma(200,shape=2,scale=2),rgamma(200,shape= scale=混合分布密dfn=function(x,a,alpha1,alpha2,theta) a*dgamma(x,shape=alpha1,scale=theta)+(1-a) dgamma(x,shape=alpha2,scale=+混合分布函pfn=function(x,a,alpha1,alpha2,theta) a*pgamma(x,shape=alpha1,scale=theta)+(1-a) pgamma(x,shape=alpha2,scale=+使用mde估計(jì)混合分布的參數(shù)，對(duì)于組數(shù)據(jù)只能使用CvM法mde.est4=mde(dat,pfn,start=list(a=0.4,alpha1=1,alpha2= theta=2.5),measure=(para= alpha1alpha2 2.0552 plot(density(dat),ylim=c(0,0.1),main="fittedcurve(dfn(x,a=para[1],alpha1=para[2],alpha2= theta=para[4]),from=-8,to=40,col="red",add=legend(20,0.08,legend=c("kerneldensity","CvM"),col= "red"),lty= N=400Bandwidth=感簡(jiǎn)單起見(jiàn)，先從單參數(shù)擬合問(wèn)題開(kāi)始，這是一個(gè)一維優(yōu)化問(wèn)題。首先生成50組來(lái)自于rate=110。然后，對(duì)于每一組隨機(jī)數(shù)，分別用基于距離的估計(jì)方法和極大似然估計(jì)進(jìn)行參數(shù)估計(jì)，將50次模擬結(jié)果的均值和標(biāo)準(zhǔn)差記錄下來(lái)。之后，將1020，30…200。不斷增大樣本量，并重復(fù)之前的過(guò)程。最終得到的結(jié)———●——●●●—————●●———●——●●●—————●●●●●●●●●●●●●●— ——————————————————————————●●●●————————●●—●●●●————————●—●———————parameterparameter0.81.01.2 指數(shù)分布隨機(jī)數(shù)，再混入兩個(gè)來(lái)自0,300]均勻分布的隨機(jī)數(shù)，再重新對(duì)參數(shù)進(jìn)行估計(jì)，結(jié)————————— ———●————————— ———●——●●●●●●● ●●●●———————— ————————●——●●●●●—●●●●parameterparameter 對(duì)于兩參數(shù)的估計(jì)，由于某些分布要求參數(shù)恒為正，使用mde函數(shù)經(jīng)常會(huì)報(bào)錯(cuò)，通常的解θ=exp(τ τ(?∞∞)θτ的值后，再例子pgammalog=function(x,logshape,logscale) pgamma(x,exp(logshape),+aa=rgamma(200,shape=3,scale=estlog=mde(aa,pgammalog,start=list(logshape=1.3,logscale= measure="CvM",method="L-BFGS-B",lower=c(0.5,- upper=c(1,5, 損失隨量的修我們知道，由于某些保險(xiǎn)條款規(guī)則的存在，保險(xiǎn)實(shí)際賠付額往往和實(shí)際損失數(shù)額是不相等的。隨量X為Y么Y=fX，f包免賠(deductible)：損失額超過(guò)某個(gè)數(shù)額，則賠付超出的部分。具體可以分為一般免陪額(ordinarydeductible)和絕對(duì)免賠(franchisedeductible)。一般免賠額的數(shù)學(xué)形式Y(jié)=(X?+

X≤dX? X>

絕對(duì)免賠額的數(shù)學(xué)形式 X≤Y X>

最大保障損失(umcoveredloss)：是保險(xiǎn)人對(duì)于單個(gè)損失支付的最大賠付額。無(wú)論數(shù)學(xué)形Y=X∧u

X≤ X>

其中隨量X∧u稱(chēng)為有限損失隨量，其期望E(X∧u)稱(chēng)為有限期望值(也就是前l(fā)ev)u(policylimit)與最大保障損失的區(qū)別在于：當(dāng)存在免賠額時(shí)，保單限額=最大保障損失-免賠額。因此當(dāng)不通貨膨脹和共同保險(xiǎn)(coinsurance)：通貨膨脹是指未來(lái)的賠付額等于當(dāng)前損失額乘以一個(gè)通脹因子，Y=1+rX；而共同保險(xiǎn)是指對(duì)每一次損失，保險(xiǎn)公司只賠付一定的比例Y=αX0<α<Y而共同保險(xiǎn)的數(shù)學(xué)形式為

(1+r)X≤(1+r)X?dd<(1+r)X≤uu?d (1+r)X>u X≤Y= α(X?d)d<X≤uα(u?d) u<X

可以看出，通貨膨脹首先對(duì)隨量X進(jìn)行通脹修正，再進(jìn)行免賠額和賠償限額的修正；共同保險(xiǎn)首先對(duì)隨量X進(jìn)行免賠額和賠償限額的修正，再進(jìn)行共保修正，兩者的順序是不YcostperlossYL)和cotperpetYP)別在于，YL是指每次損失帶來(lái)保險(xiǎn)公司的賠付額，而YP是指每次賠付帶來(lái)保險(xiǎn)公司的賠付LYPXP=YLYL=0而YPYL=(X?+

X≤dX? X>

YP

無(wú)定義X≤dX? X>

L和YPdd據(jù)要么不完全，要么干脆無(wú)法獲得，因此精算師索性將這部分損失數(shù)據(jù)忽略，只考慮導(dǎo)致正的賠付額的。如果只正的賠付么YP就個(gè)條件隨量，也以X>d為條YP=YL|X> 這樣YP的分布就是一個(gè)條件分布，而且YP的取值于0當(dāng)精算師需要根據(jù)賠付數(shù)據(jù)Y對(duì)損失額X隨隨數(shù)間XY的修正分布形式。在這個(gè)過(guò)程中，原將實(shí)際賠付數(shù)據(jù)帶入修正分布中，使用極大似然或其他估計(jì)方法估計(jì)得到修正分布的參數(shù)。根據(jù)原始分布和修正分布的關(guān)系得到原始分布的參actuar包中，coverage這個(gè)函數(shù)可以完成將原始分布變換為修正分布的工作。coverage輸出的函數(shù)語(yǔ)法coverage(pdf,cdf,deductible=0,franchise=FALSE,limit=Inf,coinsurance=1,inflation=0,per.loss=FALSE)使用說(shuō)明如果pdf和cdf同時(shí)指定，那么輸出是修正后的pdf，如果只指定cdf，那么輸出的是修正后的cdf。特別注意的是如果存在deductible或limit，那么cdf必須指定。limit設(shè)置最大保障損失u，默認(rèn)為無(wú)上限coinsurance是共保因子α，取值為0-1之間的數(shù)in?ationr0-1per.lossYPYLYPcoverage返回的是一個(gè)函數(shù)對(duì)象，如果存在重概率點(diǎn)(probabilitymass)，那么這個(gè)函數(shù)pdf并不是完全連續(xù)，probabilitymass點(diǎn)，其“分布密度”的取值其實(shí)是一個(gè)概率值，而其他點(diǎn)的取值則是密度值，在繪制對(duì)應(yīng)的分布密度圖像時(shí)，應(yīng)該對(duì)probability進(jìn)行強(qiáng)調(diào)。例子假設(shè)原始損失服從形狀參數(shù)shape=3，尺度參數(shù)scale=1的Gamma分布。首先計(jì)算修正后的密度函數(shù)，然后分別作出YP和YL的分布密度函數(shù)和分布函數(shù)。在作圖時(shí)，有以下兩點(diǎn)需 在出現(xiàn)重概率點(diǎn)時(shí)，對(duì)應(yīng)的分布函數(shù)存在跳躍par(mfrow=c(2,YP=1=7YPf=coverage(pdf=dgamma,cdf=pgamma,deductible=1,limit=curve(dgamma(x,3,1),xlim=c(0,10),ylim=c(0,0.3),ylab= main="pdfofpercurve(f(x,3,1),xlim=c(0.01,5.99),col=4,add=points(6,f(6,3,1),pch=21,bg=YL=1=7YLf1=coverage(pdf=dgamma,cdf=pgamma,deductible=1,limit= per.loss=curve(dgamma(x,3,1),xlim=c(0,10),ylim=c(0,0.3),ylab= main="pdfofpercurve(f1(x,3,1),xlim=c(0.01,5.99),col=4,add=points(6,f1(6,3,1),pch=21,bg=points(0,f1(0,3,1),pch=21,bg=YPYLF=coverage(cdf=pgamma,deductible=1,limit=curve(pgamma(x,3,1),xlim=c(0,10),ylim=c(0,1),ylab= main="cdfofpercurve(F(x,3,1),xlim=c(0,5.99),col=4,add=curve(F(x,3,1),xlim=c(6,10),col=4,add=F1=coverage(cdf=pgamma,deductible=1,limit=7,per.loss=curve(pgamma(x,3,1),xlim=c(0,10),ylim=c(0,1),ylab= main="cdfofpercurve(F1(x,3,1),xlim=c(0,5.99),col=4,add=curve(F1(x,3,1),xlim=c(6,10),col=4,add=5 cdfofper cdfofper02 802 8xxChapter風(fēng)險(xiǎn)復(fù)合分每年發(fā)生的故次數(shù)N額X從SS=X1+X2+...+ 可以看出S是一個(gè)隨機(jī)和，我們把事故次數(shù)N的分布稱(chēng)作索賠頻率分布(frequencydistribution)X(severitydistribution)S的分布稱(chēng)為復(fù)合分布(compounddistribution)?？赡芨P(guān)心這個(gè)分布。對(duì)于S的分布，我們有： FS(x)=P(S≤x)

P(S≤x|N=n)pn

XF X X其中，pn=P(N=n)是頻率分布，F(xiàn)X(x)是強(qiáng)度分布，F(xiàn)?n(x)是強(qiáng)度分布的n重卷積。如果隨量X僅在0，1，2…取值，那么n重卷積的計(jì)算方法如下：XXF?n(x)X

n=FX n=

X∑xX

F?(n?1)(x?y)fX(y)n=2,3,..連續(xù)分布的離散S4.2X的n額度通常是貨幣單位的整數(shù)倍這,依賴(lài)于研究的精確程度。定義F)為連續(xù)分布函數(shù)，x為離散化后的概率函數(shù)。目前，actuar包中的discretize函數(shù)支持四種離散化方法F(x)fx=F(x+h)?F x=a,ah,b?hcdfcdfF(x)F x=fxF(x)?F(x? x=a+h,...,離散化后的cdf總是在原cdf中點(diǎn)離散化F(a+ x=fF(x+h/2)?F(x? x=a+h,...,b?cdf正好從中間穿過(guò)離散化后的cdf無(wú)偏離散化，或者說(shuō)是局部一階矩匹配

hE(X∧a)?E(X∧a+h)+1?F x=hhfx a<x< hhE(X∧b)?E[X∧b?h]?1+F x=h離散后的分布和原分布在區(qū)間[a,b]內(nèi)有相同的取值概率和期望discretizefx函數(shù)語(yǔ)法discretize(cdf,from,to,step=1,method=c("upper","lower","rounding","unbiased"),lev,by=step,xlim=NULL)使用說(shuō)明cdfxfromtoab，也就是分布主體的范圍，stephlev只在method=“unbiased”時(shí)才指byxlim與前面參數(shù)等價(jià)，詳見(jiàn)幫助文檔例子Gamma(11)進(jìn)行離散化。fu，?，frfb分別對(duì)應(yīng)上端離散化，下端離散=="upper",from=0,to=+step==="lower",from=0,to=+step==="rounding",from=+to=5,step=fb=discretize(pgamma(x,1),method="unbiased",lev= 1),from=0,to=5,step=作出離散化后的函數(shù)圖像。函數(shù)stepfun返回一個(gè)階梯函數(shù)，di?nv是差分的逆，具體使用方法參見(jiàn)。curve(pgamma(x,1),xlim=c(0,x=seq(0,5,plot(stepfun(head(x,-1),diffinv(fu)),pch=19,,col= add=plot(stepfun(x,diffinv(fl)),pch=19,,col="red",add=plot(stepfun(head(x,-1),diffinv(fr)),pch=19,,col= add=plot(stepfun(x,diffinv(fb)),pch=19,,col="yellow",add=legend(3,0.4,legend=c("Upper","Lower","Midpoint",pgamma(x,0 x復(fù)合分布的計(jì)當(dāng)我們對(duì)索賠強(qiáng)度分布離散化后，我們可以使用一些算法計(jì)算復(fù)合分布，當(dāng)然有些方法并kaggrstPanjer(a,b,0)(a,b,1)分布族1，索賠強(qiáng)度分布正態(tài)近似法。給定頻率分布和強(qiáng)度分布就可以利用μS=E(S)=E(N)? Sσ2=Var(S)=E(N)Var(X)+Var(N S計(jì)算復(fù)合分布的均值和方差，再利F(x)≈Φ(x?μS 計(jì)算復(fù)合分布的分布函數(shù)。這種方法僅利用了前二階矩的信息，對(duì)于大樣本近似效果較正態(tài)冪近似法(NormalPower√ FS(x)=Φ(?3 9+1+6x?μS γS γS其中γS是偏度系數(shù)。正態(tài)冪近似法的原理是對(duì)標(biāo)準(zhǔn)化后的隨量S展開(kāi)為標(biāo)準(zhǔn)正態(tài)隨量及其2次冪的線性組合，也就是令S≈g(Y)，Y服從標(biāo)準(zhǔn)正態(tài)分布。這種近似法在x>μS時(shí)可以進(jìn)行，在γS<1時(shí)效果較好。因此可以對(duì)分布的右尾進(jìn)行近似。NXSSFS(x)。simul(后面會(huì)詳細(xì)講)對(duì)頻率分布和強(qiáng)度分布進(jìn)行模擬，適用于復(fù)雜系aait的參數(shù)依據(jù)所選方法的不同而不同。從下面函數(shù)的語(yǔ)法可以看出，該函數(shù)是函數(shù)語(yǔ)法"simulation"),model.freq=NULL,model.sev=NULL,p0=NULL,x.scale=1,moments,nb.simul,...,tol=1e-06,maxit=500,echo=FALSE)使用說(shuō)明遞推法中：method=“recursive”，model.freq必須是”binomial”,”geometric”,”negativebinomial”,”poisson”中的一種，分布參數(shù)可以在省略號(hào)處指定，但是參數(shù)名稱(chēng)必須要與這四個(gè)分布規(guī)定的一致。model.sevX0，1，2...個(gè)貨X0X不能取到某個(gè)值(比如2)，那么向量的對(duì)應(yīng)位置(這里是第三個(gè)元素)一定要寫(xiě)成0，通常model.sev可discretize的結(jié)果。P0N=0(ab0)分布族零調(diào)整的概率。x.scale指定X的貨幣單位，比如1元，100元等。1(a,b,0)分布族和(a,b,1)分布族是一類(lèi)滿(mǎn)足遞推關(guān)系分布的總稱(chēng)，包括泊松，二項(xiàng)，幾何和負(fù)二項(xiàng)分布及其在點(diǎn)概率調(diào)整后的分布，對(duì)該分布族的具體說(shuō)明可以參閱[5]卷積法中：method=”convolution”，model.freqN0，1，2N0的概率。model.sev的使用方法與遞推法相同。x.scaleX的貨幣單位，比1元，100元等。S正態(tài)近似法和正態(tài)冪近似法：method=”normal”或”npower”，moments是一個(gè)向量，分S的均值、方差和偏度系數(shù)，moments=c(μSσ2γS)，其中正態(tài)近似只需指定前x>μSγS<1，否則會(huì)發(fā)生報(bào)錯(cuò)。S模擬法中：method=”simulation”，具體使用方法可以參考后文中 函數(shù)的介紹nb.simul是模擬的次數(shù)函數(shù)返回的是一個(gè)aggregateDist對(duì)象，可以對(duì)其進(jìn)行五數(shù)總括(summary)，輸出結(jié)果(print)，求均值(mean)，求分位數(shù)(le)，做圖(plot),求節(jié)點(diǎn)(konts)等操作。例子X(jué)010N08，貨幣單位數(shù)設(shè)25。通過(guò)作圖觀察S的最大值是不是10×25×8=2000？par(mfrow=c(2,fx1=c(0,0.15,0.2,0.25,0.125,0.075,0.05,0.05,0.05, pn1=c(0.05,0.1,0.15,0.2,0.25,0.15,0.06,0.03,Fs1=aggregateDist("convolution",model.freq=pn1,model.sev= x.scale=Gammapoisson分布特別指定poisson分布參數(shù)lambda=10。fx2=discretize(pgamma(x,2,1),from=0,to=22,step= method="unbiased",lev=levgamma(x,2,Fs2=aggregateDist("recursive",model.freq="poisson",model.sev= lambda=10,x.scale=正態(tài)近似和正態(tài)冪近似，注意正態(tài)冪近似的有效范Fs3=aggregateDist("normal",moments=c(200,Fs4=aggregateDist("npower",moments=c(200,200,模擬法model.freq=expression(data=model.sev=expression(data=rgamma(100,Fs5=aggregateDist("simulation",nb.simul=1000, AggregateClaimAmountEmpiricalMin.1stQu. Mean3rdQu. AggregateClaimAmountExact

AggregateClaimAmountRecursive AggregateClaimAmount

AggregateClaimAmountNormal xxAggregateClaimAmountApproximation xVaR和 是常見(jiàn)的風(fēng)險(xiǎn)度量指標(biāo)，是衡量金融機(jī)構(gòu)償付能力的重要，其定義分別為VaRp=inf{x:P(S>x)≤ TVaRp=E(S|S> VaRTVaR用來(lái)計(jì)算這兩個(gè)指標(biāo)，TVaRCTE(ConditionalTailExpec-函數(shù)語(yǔ)法VaR(x,conf.level=c(0.9,0.95,0.99),names=TRUE,TVaR(x,conf.level=c(0.9,0.95,0.99),names=TRUE,使用說(shuō)明x目前僅支持aggregateDist對(duì)象，依據(jù)計(jì)算復(fù)合分布時(shí)采用的不同方法，計(jì)算VaRconf.level指定置信水平，默認(rèn)值是0.9，0.950.99names控制輸出是否包含名稱(chēng)，默認(rèn)是TURE例子考慮上文模擬法得到的結(jié)果Fs5，首先計(jì)算其VaRVaR(Fs5,names=[1]5668590198%TVaRTVaR(Fs5,conf.level=Chapter保單組復(fù)合層次模型的分布和索賠強(qiáng)度分布?；氐缴衔哪M復(fù)合分布的例子，如果認(rèn)為某保單組2010年賠付的次數(shù)服從100Poisson分布，每次賠付的賠付額服從參數(shù)(1002Gamma分布。在一次模擬中，可以首先生成一PoissonN0，然后N0Gamma分布隨機(jī)數(shù)，加總起來(lái)，就得到一個(gè)總賠付額的模擬值。如1000次，就可以得到總賠付額的經(jīng)驗(yàn)分布。()(class)，同一類(lèi)別中的不同保單tract)，以及同一保單不同的事故年度ear)，這些外生變量通過(guò)決定分布的參數(shù)進(jìn)而影響分布。因此索賠頻率和索賠強(qiáng)度的模擬可以通過(guò)復(fù)合層次模型(compundhierahicalodel)例如，考慮如下三層復(fù)合層次模型Sijt=Xijt1+···+XijtN

其中i1Ij1Jit1nij，同時(shí)jtj,Φ～Pjtj)Λji～Gamma(Φi, Φi～juj,Ψi～Lognormal(Θij,Θj|Ψi～N(Ψi, Ψi～N(2,隨量ij，Ψi和Θj在精算文獻(xiàn)中通常被稱(chēng)作風(fēng)險(xiǎn)參數(shù)；wijt是先驗(yàn)的權(quán)重。Nijtii確定后，通過(guò)參數(shù)的傳遞，可以確定jiΦiΛjΦi都確定后，就可以得到第三層分布也就Nijt的分布。jt(eposurease)都生事故1規(guī)模越大，其平均發(fā)生的損失次數(shù)也越多，如果某年車(chē)隊(duì)的車(chē)輛數(shù)為100，其風(fēng)險(xiǎn)單位數(shù)就100車(chē)年。在模型(5.2)中，索賠頻率分布服從Poisson分布，其參數(shù)的形式是jtΛ，wijt反映了在特定的風(fēng)險(xiǎn)類(lèi)別下(ij)(t)風(fēng)險(xiǎn)單位數(shù)的大小Λj是不同風(fēng)險(xiǎn)類(lèi)別(ij)所導(dǎo)致的對(duì)基本風(fēng)險(xiǎn)單位的調(diào)整因子。風(fēng)險(xiǎn)單位數(shù)越大，調(diào)整因子越大，那么對(duì)應(yīng)的Poisson分布的均IJt都是已經(jīng)確定好的數(shù)值，比如I=2J1=4J2=3n11=...=n14=4n21=n22=n23=5，如果我們把IJ，t分別與前面的風(fēng)險(xiǎn)類(lèi)別、保單和事故年度相對(duì)應(yīng)，那么上面的賦值就可以解讀為(參見(jiàn)下面的樹(shù)形示意圖)：我們要模擬一個(gè)保單組214231中的每份保單保障時(shí)間都為4年，類(lèi)別2中的每份保單的保障年度都為5年，因此我們總共443531次索賠頻(每份保單每一年都要進(jìn)行模擬)，而且對(duì)于每一次索賠，還要分叉處我們稱(chēng)之為節(jié)(nodes)，也就是對(duì)應(yīng)IJt，所謂層次模型就是通過(guò)對(duì)參數(shù)的層層遞歸，從逐步將隨機(jī)生成的參數(shù)傳遞到低層，來(lái)反映風(fēng)險(xiǎn)類(lèi)別、保單和事故年度的影響。函數(shù)simpf(simulationportfolio，也可以寫(xiě)作simul，兩個(gè)函數(shù)等價(jià))可以對(duì)復(fù)合層次模型模擬。關(guān)于復(fù)合層次模型及其模擬的更詳細(xì)的介紹，請(qǐng)參閱文獻(xiàn)[4]。函數(shù)語(yǔ)法simul(nodes,model.freq=NULL,model.sev=NULL,weights=或simpf(nodes,model.freq=NULL,model.sev=NULL,weights=使用說(shuō)明nodeslistlist對(duì)象的每一個(gè)元素是一個(gè)數(shù)值向量，每個(gè)數(shù)值向量指定在該層次下的節(jié)點(diǎn)數(shù)，需要按由到低層的順序指定，比如I，J，t的順序。model.freqmodel.sevR(expression)R(rgamma)，weightsIJt例子I，J，tnodes=list(class=2,contract=c(4,3),year=c(4,4, 4,5,5,1(I1的第一份(J1)的各(n11=4)是一的向量，緊1(I1(J2的各(n12=4)，也是一個(gè)4的向量，依此類(lèi)推，每張保單每一年有一個(gè)權(quán)重，因此共需要31個(gè)。簡(jiǎn)便起見(jiàn)，權(quán)重采用0.52.5wijt=runif(31,0.5,索賠頻率分布的模型。其中class對(duì)應(yīng)i，ct對(duì)應(yīng)Λji，r對(duì)應(yīng)jtj,i。我們可以對(duì)參數(shù)進(jìn)行更為靈活的指定，比如對(duì)參數(shù)進(jìn)行數(shù)學(xué)變換contract=rgamma(log(class),year=rpois(weights*contract*class)mf=expression(class=rexp(2),contract=rgamma(class, year=rpois(weights*索賠強(qiáng)度分布的模型ms=expression(class=rnorm(2,sqrt(0.1)),contract= 1),year=rlnorm(contract,對(duì)構(gòu)建好的復(fù)合層次模型進(jìn)行實(shí)際模擬pf=simpf(nodes=nodes,model.freq=mf,model.sev=ms,weights=PortfolioofclaimFrequencymodel ~rexp(2)contract~rgamma(class, ~rpois(weights*Severity ~rnorm(2,contract~rnorm(class,1) ~rlnorm(contract,1)Numberofclaimsper110011120230131101144341211210422151352300220函數(shù)simpf將模擬好的結(jié)果對(duì)象在pf中，默認(rèn)輸出索賠頻率和索賠強(qiáng)度的模型，以及一份保單在某一年的發(fā)生的索賠數(shù)量。比如，上面結(jié)果中第二行第四列的數(shù)字表示在類(lèi)別1的第二份保單，其第二年發(fā)生的數(shù)為2。模擬結(jié)果的處在損失數(shù)據(jù)模擬好之后，我們可以通過(guò)函aggregate,frequency，severityweights來(lái)調(diào)函數(shù)語(yǔ)法aggregate(x,by=names(x$nodes),FUN=sum,classification=TRUE,prefix=NULL,...)severity(x,by=head(names(x$node),-1),splitcol=NULL,classification=TRUE,prefix=NULL,...)weights(object,classification=TRUE,prefix=NULL,使用說(shuō)明x是一個(gè)portfolio對(duì)象，通常是simul的返回值byclass，contractyear等，可以一次指FUN(sum)(length)(mean)，中位數(shù)(median)，最大值(max)，最小值(min)等。pre?x為被匯總列添加前 提取某一年的賠付次數(shù)，具體使用方法見(jiàn)例例子aggregate返回某保單在某一年的總索賠數(shù)額，也就是Sijt。的原理類(lèi)似于Excel中的數(shù)據(jù)表11121314212223按保單類(lèi)別和保單年度分別進(jìn)行匯總，返回某類(lèi)別保單在某一年的平均索賠額度。比如，第一個(gè)類(lèi)別總共有四張保單，其中第三張保單在第一年發(fā)生一次索賠，總的賠付額為5.5，第四張保單在第一年發(fā)生四次索賠，總的賠付額為.90。因此第一類(lèi)別第一年度平均賠付=(9.908+5.251)/(1+4)=3.032aggregate(pf,by=c("class","year"),FUN= year.2year.3year.4year.5 2301.317123.50488.289 函數(shù)frequency返回某份保單在某一年的發(fā)生的索賠數(shù)量，也就是simpf默認(rèn)輸出的第aggregateaggregateFUN=length時(shí)，兩者是等價(jià)的。同樣，我們可以通過(guò)by參數(shù)進(jìn)行分門(mén)別類(lèi)的匯總。frequency(pf,prefix=110011120230131101144341211210422151352300220frequency(pf,by=12最后，函數(shù)eeri(f)返回每次賠付的索賠數(shù)額，也就是式(51)中的Xju。由于某個(gè)類(lèi) 次數(shù)小于最大次數(shù)時(shí)，空缺的數(shù)據(jù)用NA表示。11121314212223splitcol來(lái)提取出發(fā)生在某一年的賠付。比如我們要把第一個(gè)事故年度發(fā)生的賠案提取出來(lái)，severity(pf,splitcol=1112131421222311121314212223把前兩個(gè)事故年度發(fā)生的賠案提取出來(lái)，在$split處，剩余結(jié)果保存$mainseverity(pf,splitcol=c(1,classcontractclaim.1claim.2claim.3claim.4claim.5claim.61112131421222 11121314212223函數(shù)weights返回每份保單在每個(gè)年度的權(quán)重，通過(guò)和wijt進(jìn)行對(duì)比可以更好的理解權(quán)賦值的順序11121314212223aggregate(pf,classification=FALSE,)/weights(pf,classification=Chapter信度信度理論簡(jiǎn)干年的損失經(jīng)驗(yàn)，那么根據(jù)該損失經(jīng)驗(yàn)是否高于或低于手冊(cè)費(fèi)率M可以適當(dāng)增加或降低費(fèi)率，以π=Z·ˉ+(1?Z)· 其中πZˉ是過(guò)去損失經(jīng)驗(yàn)的平均，M是手冊(cè)費(fèi)率。信度因子Z反映了過(guò)去損失數(shù)據(jù)的程度，一方面我們要考慮數(shù)據(jù)量的大小，另一方面要考慮損失nX1,X2,...,Xnπ就是要n1ˉXnn1期保費(fèi)的無(wú)偏估計(jì)，而信度理論告訴我們最優(yōu)的估計(jì)是對(duì)樣本均值和手冊(cè)費(fèi)率的一個(gè)，是第n+1n信度通?？梢苑譃橛邢薏▌?dòng)信度和最精確信度。有限波動(dòng)信度建立在傳統(tǒng)的統(tǒng)計(jì)學(xué)的內(nèi)，類(lèi)似于抽樣理論中在一定的相對(duì)誤差和置信水平下確定樣本量的思想，通過(guò)對(duì)

正態(tài)似，得到“完全”條件下所需要的樣本量，也就是經(jīng)驗(yàn)損失ˉ完全可以作為下一期保費(fèi)0 條件時(shí)，通過(guò)平方根準(zhǔn)則Z=√n/n確定信度因子。最精確信度包0B¨hlmannB¨hlmannStraub模型，通過(guò)構(gòu)造過(guò)去經(jīng)驗(yàn)損失X1X2Xn的線性組合，n1[5]和[2層次信度模型(HierarhicalCreiiliyModel)由eell在1975年節(jié)自于不同風(fēng)險(xiǎn)類(lèi)別的并行數(shù)據(jù)，只要這些損失數(shù)據(jù)具有某種層次結(jié)構(gòu)，4.1小節(jié)就是講述如何模擬這種損失數(shù)據(jù)的。事實(shí)上，保險(xiǎn)公司可以用整個(gè)業(yè)務(wù)類(lèi)別(比如火險(xiǎn))的損失數(shù)據(jù)，只要這種考慮一個(gè)兩層的信度模型：在一個(gè)保單組合中，損失數(shù)據(jù)按照不同的類(lèi)(class)進(jìn)行劃分，每個(gè)類(lèi)別包含不同的合(contract)Xijt為每風(fēng)險(xiǎn)單位的損(該年度總損Sijt/wijt)，腳標(biāo)i1I表示類(lèi)別j=1Ji表示i個(gè)類(lèi)別的不同合同t=1nij表示不同年度的觀測(cè)。每個(gè)觀測(cè)都對(duì)應(yīng)一個(gè)權(quán)重wijt，表示對(duì)應(yīng)的風(fēng)險(xiǎn)單位數(shù)。在4.1節(jié)，我們介紹了如何模擬具有這種層次結(jié)構(gòu)的數(shù)據(jù)，不過(guò)在那一小節(jié)，我們稱(chēng)之為三層(各年t的索賠頻率和強(qiáng)度作為一層)。我們分別用隨量Φi和Θj代表不同類(lèi)別和合同的風(fēng)險(xiǎn)水平，模型有如下假設(shè)：隨量1,...,ΦI獨(dú)立同分布隨量1,...,ΘiJi條件獨(dú)立，給定ij對(duì)于所有tu=1Ejtj,Φi]=μ(Θij,σ2ji) t=Cov(Xijt,Xiju) t?=定義如下保單組合的結(jié)構(gòu)參數(shù)，首先是聚合保費(fèi)μ=E[E[μ(Θij,ii 合同內(nèi)方差的期望合同間方差(類(lèi)別內(nèi)方差

v=E[E[σ2(Θij, a=E[Var[μ(Θij,ii 類(lèi)別間方差

b=Var[E[μ(Θij,ii] ，B¨hlmannB¨hlmann們的目標(biāo)就是根據(jù)歷史的損失經(jīng)驗(yàn)，估計(jì)每份合同的風(fēng)險(xiǎn)保費(fèi)μ(Θij,iμ(Φi)=E[μ(Θij,ii?ij zijXijw+(1?zij)?i π?i=ziXizw+(1?zi)? 其中的信度因子為

j

wijΣ+

ziΣ+

數(shù)據(jù)的平均為X

=∑wijtX

=

zij

t=1

μvabμI ?=Xzzw iXizw zΣ

Ji

1

∑∑

wijt(Xijt?Xijw I

(nij—

i=1j=1有三種估計(jì)a和b的方法，關(guān)于這三種估計(jì)方法的更詳細(xì)請(qǐng)參見(jiàn)迭代偽估計(jì)量(iterativepseudoestimator) ?= ∑∑z

— i=1(Ji—

i=1

?b

1

zi(Xizw—Xzzw

I?

i=1這些估計(jì)量看起來(lái)簡(jiǎn)單直觀。之所以a??b稱(chēng)為偽估計(jì)量是因?yàn)榈仁接覀?cè)的信度因子需要假設(shè)μ,vab已知，而在估計(jì)時(shí)這些量只能使用其估計(jì)值。方程兩邊的相互決定的特性也決定了只能使B¨hlmannGisler估計(jì)(B¨hlmannGislerestimator)。首先記 ∑Ai wijΣ(Xijw—Xiww)2?(Ji— ci=wiΣΣ

zB ˉ d=zΣΣ

∑

∑ zΣ zz

i=1E(AiciaE(Bdb，因此有B¨hlmannGisler估計(jì)

1I

max(ci, ,?b=max(B ,dB¨hlmannGisler估計(jì)量首先在0處截?cái)?，然后再取平均，因此是有偏的估?jì)Ohlsson估計(jì)(Ohlssonestimator) II

i=1d

與B¨hlmann–Gisler原理類(lèi)似，只不過(guò)Ohlsson估計(jì)量將求平均方式改為平均。a?′′簡(jiǎn)化了中間節(jié)點(diǎn)的計(jì)算形式actuarcm函數(shù)支持這三種方法擬合層次信度模型，cmcredibilitymodel的縮cmRlmcm”，該包的作者在幫助文檔和發(fā)表的文章中不遺余力地重復(fù)著這個(gè)冷笑話(huà)cmB¨hlmann，B¨hlmann函數(shù)語(yǔ)法cm(formula,data,ratios,weights,subset,regformula=NULL,regdata,ercept=FALSE,method=c("Buhlmann-Gisler","Ohlsson","iterative"),tol=sqrt(.Machine$double.eps),maxit=100,echo=FALSE)使用說(shuō)明formulalm中的公式寫(xiě)法，但是?左邊不需指定因變量。?左邊按順序指定分層的因子，比如一個(gè)保單組合分為不同的類(lèi)別(class)，在每個(gè)(contract)，如果令:+表示不同項(xiàng)的分割，那么可以令?rmula=?class+class:contract。Xijtwijt(B¨hlmann模ratiosweights指明哪些列是Xijt或wijt，weights取觀測(cè)的子集進(jìn)行模型擬合，需要輸入邏輯表達(dá)式，要使該參數(shù)有 必須是數(shù)據(jù)(data.frame)regformula，regdata和ercept是信度回歸模型的參數(shù)，將在下一小節(jié)講methodab參數(shù)的估計(jì)方法，默認(rèn)是Buhlmann-Gisler法。當(dāng)層數(shù)為1B¨hlmannB¨hlmannStraub時(shí)，三種方法估計(jì)結(jié)果相tol指定迭代收斂的條件，maxit指定最大迭代次數(shù)，echo指定是否輸出迭代過(guò)程，默認(rèn)例子B¨hlmandat1=data.frame(plyhder=c(1,2),ratio.1=c(NA, ratio.2=c(10000/50,21000/110),ratio.3=c(13000/60, weight.1=c(NA,100),weight.2=c(50,110),weight.3= fit1=cm(~plyhder,dat1,ratios=ratio.1:ratio.3,weights=cm(formula=~plyhder,data=dat1,ratios=ratio.1:ratio.3,weights=weight.1:weight.3)StructureParametersEstimatorsCollectivepremium:191.7Betweenplyhdervariance:380.9Withinplyhdervariance:17831果

cm(formula=~plyhder,data=dat1,ratios=ratio.1:ratio.3,weights=weight.1:weight.3)StructureParametersEstimatorsCollectivepremium:191.7Betweenplyhdervariance:380.9Withinplyhdervariance:17831DetailedpremiumsLevel:plyhderplyhderIndiv.meanWeightCred.factorCred. 203.9如果要擬合B¨hlmann模型，可以不指定權(quán)重，但每張保單的觀測(cè)年數(shù)要相dat2=data.frame(plyhder=c(1,2),ratio.1=c(10000/50, ratio.2=c(13000/60,21000/110),ratio.3=c(10200/52,fit2=cm(~plyhder,dat2,ratios=ratio.1:ratio.3,method=cm(formula=~plyhder,data=dat2,ratios=ratio.1:ratio.3,method="iterative")StructureParametersEstimatorsCollectivepremium:190.9Betweenplyhdervariance:300.0Withinplyhdervariance:166.8下面擬合層次信度模型，數(shù)據(jù)我們采 函數(shù)進(jìn)行模擬，索賠強(qiáng)度和索賠頻率模型4.1wijt[50100]wijt=runif(31,80,nodes=list(class=2,contract=c(4,3),year=c(4,4, 4,5,5,mf=expression(class=rexp(1),contract=rgamma(class+ 1),year=rpois(weights*ms=expression(class=rnorm(2,sqrt(0.1)),contract= 1),year=rlnorm(contract,pf=simpf(nodes=nodes,model.freq=mf,model.sev=ms,weights=計(jì)算Xijt=Sijt/wijt，并構(gòu)造數(shù)ratio.mat=aggregate(pf,classification=FALSE,prefix= dat3=cbind(weights(pf,prefix="weight."),擬合層次信度模型fit3=cm(~class+class:contract,data=dat3,ratio= weights=cm(formula=~class+class:contract,data=dat3,ratios=ratio.year.1:ratio.year.5,weights=weight.year.1:weight.year.5)StructureParametersEstimatorsCollectivepremium:23.3Betweenclassvariance:391.3Withinclass/Betweencontractvariance:213.2Withincontractvariance:630.3DetailedpremiumsLevel:classclassIndiv.meanWeightCred.factorCred.

contractIndiv.meanWeightCred.factorCred.11121314212223[1]36.29[1]52.330 7.785 5.612還可以按照不同層次進(jìn)行輸出結(jié)果和summary(fit3,levels=cm(formula=~class+class:contract,data=dat3,ratios=ratio.year.1:ratio.year.5,weights=weight.year.1:weight.year.5)StructureParametersCollectivepremium:Betweenclassvariance:391.3Withinclassvariance:213.2DetailedpremiumsLevel:classclassIndiv.meanWeightCred.factorCred. predict(fit3,levels=[1]52.330 7.785 5.612信度回歸模B¨hlmannStraub模型中，需要假設(shè)各年損失數(shù)據(jù)是平穩(wěn)的，沒(méi)有明顯趨勢(shì)，如果數(shù)據(jù)呈現(xiàn)系統(tǒng)性的趨勢(shì)，那么B¨hlmann–Straub模型就會(huì)低估或者高估真實(shí)保費(fèi)。精算師Hachemeister就遇到了這種情況，他想要各州第責(zé)任險(xiǎn)在未來(lái)一年的平均賠付額度，B¨hlmannStraub模型中的平穩(wěn)性假設(shè)Eji關(guān)于j獨(dú)立)遭到了破壞。Hachemeister在其1975年的著名文章中首先介紹一下在這篇文章中所使用的數(shù)據(jù)。該數(shù)據(jù)包含了5個(gè)州從1970年7月1973年6月共12個(gè)季度的機(jī)動(dòng)車(chē)第責(zé)任險(xiǎn)的損失數(shù)據(jù)。數(shù)據(jù)總共有5行25列，每行代一個(gè)州，第1列為每個(gè)州的，第2至13列為某個(gè)州在某季度的案均賠付額，第14至列為對(duì)應(yīng)的數(shù)量。 中可以通過(guò)如下語(yǔ)句調(diào)用此數(shù)據(jù)12345weight.10 的的。當(dāng)然現(xiàn)間的差行111111 111111 414444444loss0 0i現(xiàn)在我們正式引入信度回歸模型。給定具有I個(gè)風(fēng)險(xiǎn)的保單組合，令X′=(Xi1,...,Xin)iwi′=(wi1win)為對(duì)應(yīng)的權(quán)重向量。比如在Hachemeister數(shù)據(jù)中，Xi代表第i個(gè)州的按均賠付額，wi代表第i個(gè)州的數(shù)量。i如果令Θi代表第i個(gè)州的風(fēng)險(xiǎn)水平，因此需要在給定Θi的條件下，擬合回歸方程X=Yβ+ YΘiβ具有隨機(jī)效應(yīng)，是關(guān)于Θi的隨量，不妨記作β(Θi)。模型具有如下假設(shè)：ΘiXijE(Xi|Θi)= Var(Xij|Θi)隨機(jī)向量(Θ1X1)(ΘIXI)相互獨(dú)立，且ΘiΘ獨(dú)立同分布

根據(jù)條件獨(dú)立性，性質(zhì)(6.18)可以得到數(shù)據(jù)的協(xié)方差矩Cov(Xi,X′|Θi)= Wi=diag(wi1,..., 當(dāng)不存在協(xié)變量時(shí)，亦Xi11)′時(shí)1E(Xi|Θi) 1B¨hlmannE(Xij|Θi)=β0(Θi)+β1(Θi)· 寫(xiě)成矩陣形式，有

E(X|Θ)

)

β(Θ 根據(jù)如上假設(shè)，我們得到β(Θi)^(Θi)=AiBi+(I? 其中

Ai=T(T+i是信度矩陣Xi的最優(yōu)線性無(wú)偏估計(jì)。

Bi=(Y′V?1Yi)?1Y′i iiVi=E[Cov(Xi,X′iT=Cov(β(Θi),是回歸系數(shù)協(xié)方差陣(組間協(xié)方差矩陣)β=是聚合回歸系數(shù)(聚合保費(fèi))在函數(shù)cm中，regformula，ercept三個(gè)參數(shù)控制信度回歸模型。r