用r-也能做精算actuar包學習筆記

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

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

cj?1+2

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

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

1n

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

1∑n

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

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

min(xj, 在此略去，有的同學可以參考幫助文檔。uu計算出的經(jīng)驗有限期望值是不一樣的。elevu的經(jīng)驗有限期望值，只需要特別指定u即可。例子組數(shù)據(jù)，返回的lev是上限u的函數(shù)lev=這里將保單限額u200[1]lev函數(shù)的拐點，注意到拐點都發(fā)生在數(shù)據(jù)點，而觀察下圖可知拐點間的函數(shù)都是線 46107141259317351567分組數(shù)lev2=[1]分別對組數(shù)據(jù)和分組數(shù)據(jù)的有限期望函數(shù)作圖par(mfrow=c(1,plot(lev,type="o",pch=plot(lev2,type="o",pch=●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●Empirical200Empirical100 100 Chapter損失損失分布種根據(jù)損失額的特征，損失分布常選用具有非負支集(密度函數(shù)f(x)的支集指的是使得f(x)?=0x的集合)的連續(xù)分布。Rd，p，q，r四種函數(shù)，分別是密度函數(shù)、分布函數(shù)、分布函數(shù)的反函數(shù)(分位數(shù))和生成該分布的隨機數(shù)。actuar包提供了與[5]的附錄A中所列示的連續(xù)分布族相配套的這四種函數(shù)(除去逆和對數(shù)t分布，但包括Gamma分布)Rstats中并不自帶，但有些分布在精算研究中卻很(pareto分布)。此外，actuarm、lev和mgf三種函數(shù)，m是計算理論原點矩，lev是計算有限期望值，mgf是計算矩母函數(shù)。密度函數(shù)、分布函數(shù)、原點矩、有限期望值及其k次方都可以通過查詢該附錄得到。對于經(jīng)驗數(shù)據(jù)，如上面所介紹的，actuaremmelev來計算經(jīng)驗原點矩和經(jīng)驗有限期望值(這兩個函數(shù)的前綴都是empirical)。ratescale參數(shù)，scale=1/rate，因此兩者在本質上是等價的，[5]中使用的是scale參數(shù)，在指定參數(shù)時千萬不要弄混。例子這里以雙參數(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ù)可以進行極大似然估計。在actuar包中，mde函數(shù)則提供了三種距離最小化的分布擬合方法(miniumdistanceestimates)。Cram′er-vonMises(CvM)(對于分組數(shù)據(jù)是ogive)的距離。未分組數(shù)據(jù)d(θ)

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

d(θ)

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

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

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

wj(LAS(cj?1,cj;θ)?LA?Sn(cj?1, L?Sn(x,?n?njjstatsoptim函數(shù)語法mde(x,fun,start,measure=c("CvM","chi-使用說明x是分組數(shù)據(jù)對象的或未分組的數(shù)據(jù)fun是待擬合的分布，CvM法和修正卡方法需要指定分布函數(shù)：p**。LAS法需要指定理論有限期望函數(shù)lev**。start指定參數(shù)初始值。形式必須以列表的形式，形式可以見例子，有幾個參數(shù)就要指定measure是指定方法。weight指定權重，否則采用默認權optimL-BFGS-B方法進行優(yōu)化可以添加參數(shù)method=“L-BFGS-B”。mde(list)，rate是參數(shù)估計結果，distance是最小化后的距離。我們可以對上面的gdental數(shù)據(jù)進行分布擬合，這是一個分組數(shù)據(jù)，在為參數(shù)估計賦初始值時，假設數(shù)據(jù)來自均值為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= 我們還可以對組數(shù)據(jù)進行分布擬合，下面的一個例子是一個混合分布例子首先生成400個隨機數(shù)，其中200個來自Gamma(α2θ2)，200個來自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估計混合分布的參數(shù)，對于組數(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=感簡單起見，先從單參數(shù)擬合問題開始，這是一個一維優(yōu)化問題。首先生成50組來自于rate=110。然后，對于每一組隨機數(shù)，分別用基于距離的估計方法和極大似然估計進行參數(shù)估計，將50次模擬結果的均值和標準差記錄下來。之后，將1020，30…200。不斷增大樣本量，并重復之前的過程。最終得到的結———●——●●●—————●●———●——●●●—————●●●●●●●●●●●●●●— ——————————————————————————●●●●————————●●—●●●●————————●—●———————parameterparameter0.81.01.2 指數(shù)分布隨機數(shù)，再混入兩個來自0,300]均勻分布的隨機數(shù)，再重新對參數(shù)進行估計，結————————— ———●————————— ———●——●●●●●●● ●●●●———————— ————————●——●●●●●—●●●●parameterparameter 對于兩參數(shù)的估計，由于某些分布要求參數(shù)恒為正，使用mde函數(shù)經(jīng)常會報錯，通常的解θ=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, 損失隨量的修我們知道，由于某些保險條款規(guī)則的存在，保險實際賠付額往往和實際損失數(shù)額是不相等的。隨量X為Y么Y=fX，f包免賠(deductible)：損失額超過某個數(shù)額，則賠付超出的部分。具體可以分為一般免陪額(ordinarydeductible)和絕對免賠(franchisedeductible)。一般免賠額的數(shù)學形式Y=(X?+

X≤dX? X>

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

最大保障損失(umcoveredloss)：是保險人對于單個損失支付的最大賠付額。無論數(shù)學形Y=X∧u

X≤ X>

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

(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

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

X≤dX? X>

YP

無定義X≤dX? X>

L和YPdd據(jù)要么不完全，要么干脆無法獲得，因此精算師索性將這部分損失數(shù)據(jù)忽略，只考慮導致正的賠付額的。如果只正的賠付么YP就個條件隨量，也以X>d為條YP=YL|X> 這樣YP的分布就是一個條件分布，而且YP的取值于0當精算師需要根據(jù)賠付數(shù)據(jù)Y對損失額X隨隨數(shù)間XY的修正分布形式。在這個過程中，原將實際賠付數(shù)據(jù)帶入修正分布中，使用極大似然或其他估計方法估計得到修正分布的參數(shù)。根據(jù)原始分布和修正分布的關系得到原始分布的參actuar包中，coverage這個函數(shù)可以完成將原始分布變換為修正分布的工作。coverage輸出的函數(shù)語法coverage(pdf,cdf,deductible=0,franchise=FALSE,limit=Inf,coinsurance=1,inflation=0,per.loss=FALSE)使用說明如果pdf和cdf同時指定，那么輸出是修正后的pdf，如果只指定cdf，那么輸出的是修正后的cdf。特別注意的是如果存在deductible或limit，那么cdf必須指定。limit設置最大保障損失u，默認為無上限coinsurance是共保因子α，取值為0-1之間的數(shù)in?ationr0-1per.lossYPYLYPcoverage返回的是一個函數(shù)對象，如果存在重概率點(probabilitymass)，那么這個函數(shù)pdf并不是完全連續(xù)，probabilitymass點，其“分布密度”的取值其實是一個概率值，而其他點的取值則是密度值，在繪制對應的分布密度圖像時，應該對probability進行強調。例子假設原始損失服從形狀參數(shù)shape=3，尺度參數(shù)scale=1的Gamma分布。首先計算修正后的密度函數(shù)，然后分別作出YP和YL的分布密度函數(shù)和分布函數(shù)。在作圖時，有以下兩點需 在出現(xiàn)重概率點時，對應的分布函數(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ā)生的故次數(shù)N額X從SS=X1+X2+...+ 可以看出S是一個隨機和，我們把事故次數(shù)N的分布稱作索賠頻率分布(frequencydistribution)X(severitydistribution)S的分布稱為復合分布(compounddistribution)?？赡芨P心這個分布。對于S的分布，我們有： FS(x)=P(S≤x)

P(S≤x|N=n)pn

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

n=FX n=

X∑xX

F?(n?1)(x?y)fX(y)n=2,3,..連續(xù)分布的離散S4.2X的n額度通常是貨幣單位的整數(shù)倍這,依賴于研究的精確程度。定義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中點離散化F(a+ x=fF(x+h/2)?F(x? x=a+h,...,b?cdf正好從中間穿過離散化后的cdf無偏離散化，或者說是局部一階矩匹配

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]內有相同的取值概率和期望discretizefx函數(shù)語法discretize(cdf,from,to,step=1,method=c("upper","lower","rounding","unbiased"),lev,by=step,xlim=NULL)使用說明cdfxfromtoab，也就是分布主體的范圍，stephlev只在method=“unbiased”時才指byxlim與前面參數(shù)等價，詳見幫助文檔例子Gamma(11)進行離散化。fu，?，frfb分別對應上端離散化，下端離散=="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返回一個階梯函數(shù)，di?nv是差分的逆，具體使用方法參見。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復合分布的計當我們對索賠強度分布離散化后，我們可以使用一些算法計算復合分布，當然有些方法并kaggrstPanjer(a,b,0)(a,b,1)分布族1，索賠強度分布正態(tài)近似法。給定頻率分布和強度分布就可以利用μS=E(S)=E(N)? Sσ2=Var(S)=E(N)Var(X)+Var(N S計算復合分布的均值和方差，再利F(x)≈Φ(x?μS 計算復合分布的分布函數(shù)。這種方法僅利用了前二階矩的信息，對于大樣本近似效果較正態(tài)冪近似法(NormalPower√ FS(x)=Φ(?3 9+1+6x?μS γS γS其中γS是偏度系數(shù)。正態(tài)冪近似法的原理是對標準化后的隨量S展開為標準正態(tài)隨量及其2次冪的線性組合，也就是令S≈g(Y)，Y服從標準正態(tài)分布。這種近似法在x>μS時可以進行，在γS<1時效果較好。因此可以對分布的右尾進行近似。NXSSFS(x)。simul(后面會詳細講)對頻率分布和強度分布進行模擬，適用于復雜系aait的參數(shù)依據(jù)所選方法的不同而不同。從下面函數(shù)的語法可以看出，該函數(shù)是函數(shù)語法"simulation"),model.freq=NULL,model.sev=NULL,p0=NULL,x.scale=1,moments,nb.simul,...,tol=1e-06,maxit=500,echo=FALSE)使用說明遞推法中：method=“recursive”，model.freq必須是”binomial”,”geometric”,”negativebinomial”,”poisson”中的一種，分布參數(shù)可以在省略號處指定，但是參數(shù)名稱必須要與這四個分布規(guī)定的一致。model.sevX0，1，2...個貨X0X不能取到某個值(比如2)，那么向量的對應位置(這里是第三個元素)一定要寫成0，通常model.sev可discretize的結果。P0N=0(ab0)分布族零調整的概率。x.scale指定X的貨幣單位，比如1元，100元等。1(a,b,0)分布族和(a,b,1)分布族是一類滿足遞推關系分布的總稱，包括泊松，二項，幾何和負二項分布及其在點概率調整后的分布，對該分布族的具體說明可以參閱[5]卷積法中：method=”convolution”，model.freqN0，1，2N0的概率。model.sev的使用方法與遞推法相同。x.scaleX的貨幣單位，比1元，100元等。S正態(tài)近似法和正態(tài)冪近似法：method=”normal”或”npower”，moments是一個向量，分S的均值、方差和偏度系數(shù)，moments=c(μSσ2γS)，其中正態(tài)近似只需指定前x>μSγS<1，否則會發(fā)生報錯。S模擬法中：method=”simulation”，具體使用方法可以參考后文中 函數(shù)的介紹nb.simul是模擬的次數(shù)函數(shù)返回的是一個aggregateDist對象，可以對其進行五數(shù)總括(summary)，輸出結果(print)，求均值(mean)，求分位數(shù)(le)，做圖(plot),求節(jié)點(konts)等操作。例子X010N08，貨幣單位數(shù)設25。通過作圖觀察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和 是常見的風險度量指標，是衡量金融機構償付能力的重要，其定義分別為VaRp=inf{x:P(S>x)≤ TVaRp=E(S|S> VaRTVaR用來計算這兩個指標，TVaRCTE(ConditionalTailExpec-函數(shù)語法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,使用說明x目前僅支持aggregateDist對象，依據(jù)計算復合分布時采用的不同方法，計算VaRconf.level指定置信水平，默認值是0.9，0.950.99names控制輸出是否包含名稱，默認是TURE例子考慮上文模擬法得到的結果Fs5，首先計算其VaRVaR(Fs5,names=[1]5668590198%TVaRTVaR(Fs5,conf.level=Chapter保單組復合層次模型的分布和索賠強度分布。回到上文模擬復合分布的例子，如果認為某保單組2010年賠付的次數(shù)服從100Poisson分布，每次賠付的賠付額服從參數(shù)(1002Gamma分布。在一次模擬中，可以首先生成一PoissonN0，然后N0Gamma分布隨機數(shù)，加總起來，就得到一個總賠付額的模擬值。如1000次，就可以得到總賠付額的經(jīng)驗分布。()(class)，同一類別中的不同保單tract)，以及同一保單不同的事故年度ear)，這些外生變量通過決定分布的參數(shù)進而影響分布。因此索賠頻率和索賠強度的模擬可以通過復合層次模型(compundhierahicalodel)例如，考慮如下三層復合層次模型Sijt=Xijt1+···+XijtN

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

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

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

b=Var[E[μ(Θij,ii] ，B¨hlmannB¨hlmann們的目標就是根據(jù)歷史的損失經(jīng)驗，估計每份合同的風險保費μ(Θ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有三種估計a和b的方法，關于這三種估計方法的更詳細請參見迭代偽估計量(iterativepseudoestimator) ?= ∑∑z

— i=1(Ji—

i=1

?b

1

zi(Xizw—Xzzw

I?

i=1這些估計量看起來簡單直觀。之所以a??b稱為偽估計量是因為等式右側的信度因子需要假設μ,vab已知，而在估計時這些量只能使用其估計值。方程兩邊的相互決定的特性也決定了只能使B¨hlmannGisler估計(B¨hlmannGislerestimator)。首先記 ∑Ai wijΣ(Xijw—Xiww)2?(Ji— ci=wiΣΣ

zB ˉ d=zΣΣ

∑

∑ zΣ zz

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

1I

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

i=1d

與B¨hlmann–Gisler原理類似，只不過Ohlsson估計量將求平均方式改為平均。a?′′簡化了中間節(jié)點的計算形式actuarcm函數(shù)支持這三種方法擬合層次信度模型，cmcredibilitymodel的縮cmRlmcm”，該包的作者在幫助文檔和發(fā)表的文章中不遺余力地重復著這個冷笑話cmB¨hlmann，B¨hlmann函數(shù)語法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)使用說明formulalm中的公式寫法，但是?左邊不需指定因變量。?左邊按順序指定分層的因子，比如一個保單組合分為不同的類別(class)，在每個(contract)，如果令:+表示不同項的分割，那么可以令?rmula=?class+class:contract。Xijtwijt(B¨hlmann模ratiosweights指明哪些列是Xijt或wijt，weights取觀測的子集進行模型擬合，需要輸入邏輯表達式，要使該參數(shù)有 必須是數(shù)據(jù)(data.frame)regformula，regdata和ercept是信度回歸模型的參數(shù)，將在下一小節(jié)講methodab參數(shù)的估計方法，默認是Buhlmann-Gisler法。當層數(shù)為1B¨hlmannB¨hlmannStraub時，三種方法估計結果相tol指定迭代收斂的條件，maxit指定最大迭代次數(shù)，echo指定是否輸出迭代過程，默認例子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模型，可以不指定權重，但每張保單的觀測年數(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ù)進行模擬，索賠強度和索賠頻率模型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=計算Xijt=Sijt/wijt，并構造數(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還可以按照不同層次進行輸出結果和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ù)據(jù)是平穩(wěn)的，沒有明顯趨勢，如果數(shù)據(jù)呈現(xiàn)系統(tǒng)性的趨勢，那么B¨hlmann–Straub模型就會低估或者高估真實保費。精算師Hachemeister就遇到了這種情況，他想要各州第責任險在未來一年的平均賠付額度，B¨hlmannStraub模型中的平穩(wěn)性假設Eji關于j獨立)遭到了破壞。Hachemeister在其1975年的著名文章中首先介紹一下在這篇文章中所使用的數(shù)據(jù)。該數(shù)據(jù)包含了5個州從1970年7月1973年6月共12個季度的機動車第責任險的損失數(shù)據(jù)。數(shù)據(jù)總共有5行25列，每行代一個州，第1列為每個州的，第2至13列為某個州在某季度的案均賠付額，第14至列為對應的數(shù)量。 中可以通過如下語句調用此數(shù)據(jù)12345weight.10 的的。當然現(xiàn)間的差行111111 111111 414444444loss0 0i現(xiàn)在我們正式引入信度回歸模型。給定具有I個風險的保單組合，令X′=(Xi1,...,Xin)iwi′=(wi1win)為對應的權重向量。比如在Hachemeister數(shù)據(jù)中，Xi代表第i個州的按均賠付額，wi代表第i個州的數(shù)量。i如果令Θi代表第i個州的風險水平，因此需要在給定Θi的條件下，擬合回歸方程X=Yβ+ YΘiβ具有隨機效應，是關于Θi的隨量，不妨記作β(Θi)。模型具有如下假設：ΘiXijE(Xi|Θi)= Var(Xij|Θi)隨機向量(Θ1X1)(ΘIXI)相互獨立，且ΘiΘ獨立同分布

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

E(X|Θ)

)

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

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

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