« アクチュアリー試験の歴史 | トップページ | 理事長賞、大館賞、理事長特別賞 »

2012年2月19日 (日)

Rによるアクチュアリーの統計分析(第2章)

かなり許しがたい。



お断りしておきますが、アフィリエイトリンクを使っているのは書影を出すのに便利だからという理由だけで、購入をおすすめするものではありません。

先日に引き続き、「Rによるアクチュアリーの統計分析」第2章です。第2章は年金数理。

まずは年齢群団別の脱退実績データを補正し、スムーズな脱退率を作成します。とりあえず本のとおりに入力してみます(p19)。

> x<-c(22,27,32,37,42,47,52,57)
> bunbo<-c(138,67,165,154,145,135,143,152)
> bunshi<-c(9,2,6,3,2,2,2,1)
> (result1<-nls(soqw~a+b*x+c*x^2,start=c(a=0.01,b=0,c=0)))
以下にエラー nls(soqw ~ a + b * x + c * x^2, start = c(a = 0.01, b = 0, c = 0)) :
   'data' 中に初期値の無いパラメータがあります

エラーになりました。soqwが定義されていないためです。

result1の計算の前に、次の1文が抜けています。

> soqw <- bunshi/bunbo

ともかく、これを入れれば脱退率の補正のプログラムは走ります。気を取り直して先に行きましょう。

次は年金数理を勉強している方にはおなじみ、極限方程式が出てきます。定常状態の人口分布の作成例がありますが、これは正しく走ります。

定常状態ができたら、トローブリッジモデルによる財政方式の分類です。p29の計算をやってみます。

> v<-1/1.03;la<-l[1:40];lr<-l[41:88]
エラー:  オブジェクト 'l' がありません

単に本の記述をなぞっているだけなのに、なぜこんなにエラーが出るのでしょう?
ここでの変数はlではなくてeを使う必要があるようです。

> D<-(v^(0:87))*e;N<-rev(cumsum(rev(D)))
> ad<-N[41]/D[1:40];ar<-N[41:88]/D[41:88]
> aa<-(N[1:40]-N[41])/D[1:40];ib<-0:39/40
> L<-sum(la);B<-sum(lr);Sa<-sum(la*ad);Sr<-sum(lr*ar);S<-Sa+Sr
> Sap<-sum(ib*la*ad);Saf<-sum((1-ib)*la*ad);Ga<-sum(la*aa)
> Sf<-la[1]*ad[1]*sum(v*e);Gf<-la[1]*aa[1]*sum(v*e)
> c(L,B,Sa,Sr,S,Sap,Saf,Ga,Sf,Gf)
[1]    2520.8784     985.6637   16229.4639   10436.9480   26666.4119
[6]    9473.8835    6755.5804   32194.7328  732769.5312 5551444.1221

SfとGf(上の太字)が合いません。SfとGfの定義を見るに、次のように計算すべものになっています。

> Sf <- la[1]*ad[1]*sum(v^(0:87)); Gf <- la[1]*aa[1]*sum(v^(0:87))
> c(Sf, Gf)
[1]  6841.717 51832.681

ともかく、これで各財政方式の保険料を計算する要素が揃いました。さっそく計算してみます。

> (Pea<-ad[1]/aa[1])
[1] 0.1319962
> Pea*L
[1] 332.7464
> PSL<-S-Pea*Ga
> PSL/sum(vv[1:20])
エラー:  オブジェクト 'vv' がありません

またコケました。

> vv <- v^(0:87)

を追加する必要があります。

> ((Sa+Sr)/Ga)    #閉鎖型総合保険料方式
[1] 0.8282849
> ((Sa+Sr)/Ga*L)  #制度全体
[1] 2088.006
> ((Sa+Sr+Sf)/(Ga+Gf))    #開放型総合保険料方式
[1] 0.3987762
> ((Sa+Sr+Sf)/(Ga+Gf)*L)  #制度全体
[1] 1005.266

このへんは本の記述は正しいようです。

> J<-1:40
> 1/40*N[42]/D[J]  #単位積立方式保険料率
[1] 0.05043388 0.05355351 0.05686609 0.06038358 0.06411865 0.06808475
[7] 0.07229618 0.07676811 0.08151665 0.08655892 0.09191308 0.09759843
-------------------------------(省略)---------------------------------
[31] 0.25129481 0.26144814 0.27201170 0.28300207 0.29443650 0.30633292
[37] 0.31871001 0.33158718 0.34498464 0.35892342
> N[41]/(N[J]-N[41])  #個人別平準保険料率
[1]  0.1319962  0.1406197  0.1498386  0.1596983  0.1702486  0.1815433
[7]  0.1936417  0.2066083  0.2205143  0.2354375  0.2514639  0.2686883
-------------------------------(省略)---------------------------------
[31]  1.2734706  1.4450589  1.6600486  1.9370418  2.3070418  2.8258543
[37]  3.6050904  4.9051758  7.5073865 15.3181018

Rで計算した結果は上のようになったのですが、本は単位積立方式保険料率と個人別平準保険料率の数値が逆になっています。

第1章・第2章までで30ページですが、その中のプログラム記述だけで6箇所の誤植です。かつ、朝倉書店のサイトを見ても正誤表は見当たりません。

内容のレベルうんぬん以前に、ちょっとひどすぎる。

« アクチュアリー試験の歴史 | トップページ | 理事長賞、大館賞、理事長特別賞 »

アクチュアリー」カテゴリの記事

書評」カテゴリの記事

コメント

154ページのコマンドの8行目に記載の『hbe』というのは見つかりましたか?

朝倉書店のHPからコマンドとデータファイルをダウンロードしたのですが、見つかりませんでした。

大昔にお世話になった『APL』にRが非常に似ているので、懐かしさを噛みしめながら、田中先生の本と格闘しております。

> 154ページのコマンドの8行目に記載の『hbe』というのは見つかりましたか?

どこかのパッケージに入っている関数かと思ってちょっと調べてみましたが、わかりませんね。
こういう、「Rを学ぶ」のと関係ないところで苦労させられるので辟易です…

コメントを書く

コメントは記事投稿者が公開するまで表示されません。

(ウェブ上には掲載しません)

トラックバック

この記事のトラックバックURL:
http://app.cocolog-nifty.com/t/trackback/554939/54018799

この記事へのトラックバック一覧です: Rによるアクチュアリーの統計分析(第2章):

« アクチュアリー試験の歴史 | トップページ | 理事長賞、大館賞、理事長特別賞 »

フォト
無料ブログはココログ
2017年4月
            1
2 3 4 5 6 7 8
9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 29
30            

最近のトラックバック