options(width=100,warn=-1,digits=5,messages=1,tidy=TRUE)
## 倒産企業と非倒産企業の仮想データを準備
# ya1:倒産企業の自己資本比率 ya2:倒産企業のインタレスト
#カバレッジ比率

# yb1:非倒産企業の自己資本比率 yb2:非倒産企業の
#インタレストカバレッジ比率

# 便宜的に倒産企業は0 非倒産企業は1 を識別標識zに
#入力する 

ya1<-c(15.5,5.6,28.3,12.7,-78.7)
ya2<-c(1.14,1.03,-0.42,-8.05,0.58)
yb1<-c(11.4,13.7,30.25,34.06,54.24,52.23,35.61,36.47,59.7,46.8,34.3,33.4,47.7,38.1)
yb2<-c(2.3,1.93,2.43,5.5,7.76,6.05,6.31,5.2,22.36,18.92,3.21,3.53,3.91,-3.65)

#始めからデータフレーム型でデータを用意すればよいのだ
#が、線形判別分析の計算練習で、通常の多変量解析による行
#列演算で得た結果とパッケージMASSのldaを使って得た結果を#比較するために、大変に冗長な手順に
なるが,わざとベクトル#でデータを準備している。


(za1<-rep(0,length(ya1)))

## [1] 0 0 0 0 0

ya<-cbind(ya1,ya2,za1)
length(ya1)

## [1] 5

colnames(ya)<-c("x1","x2","z")
ya

##         x1    x2 z
## [1,]  15.5  1.14 0
## [2,]   5.6  1.03 0
## [3,]  28.3 -0.42 0
## [4,]  12.7 -8.05 0
## [5,] -78.7  0.58 0

(zb1<-rep(1,length(yb1)))

##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1

length(yb1)

## [1] 14yb<-cbind(yb1,yb2,zb1)
yb

##         yb1   yb2 zb1
##  [1,] 11.40  2.30   1
##  [2,] 13.70  1.93   1
##  [3,] 30.25  2.43   1
##  [4,] 34.06  5.50   1
##  [5,] 54.24  7.76   1
##  [6,] 52.23  6.05   1
##  [7,] 35.61  6.31   1
##  [8,] 36.47  5.20   1
##  [9,] 59.70 22.36   1
## [10,] 46.80 18.92   1
## [11,] 34.30  3.21   1
## [12,] 33.40  3.53   1
## [13,] 47.70  3.91   1
## [14,] 38.10 -3.65   1

colnames(yb)<-c("x1","x2","z")
yb

##          x1    x2 z
##  [1,] 11.40  2.30 1
##  [2,] 13.70  1.93 1
##  [3,] 30.25  2.43 1
##  [4,] 34.06  5.50 1
##  [5,] 54.24  7.76 1
##  [6,] 52.23  6.05 1
##  [7,] 35.61  6.31 1
##  [8,] 36.47  5.20 1
##  [9,] 59.70 22.36 1
## [10,] 46.80 18.92 1
## [11,] 34.30  3.21 1
## [12,] 33.40  3.53 1
## [13,] 47.70  3.91 1
## [14,] 38.10 -3.65 1

##以下で倒産企業と非倒産企業の説明変数全部の行列が作ら
#れる

yall<-rbind(ya,yb)
yall

##           x1    x2 z
##  [1,]  15.50  1.14 0
##  [2,]   5.60  1.03 0
##  [3,]  28.30 -0.42 0
##  [4,]  12.70 -8.05 0
##  [5,] -78.70  0.58 0
##  [6,]  11.40  2.30 1
##  [7,]  13.70  1.93 1
##  [8,]  30.25  2.43 1
##  [9,]  34.06  5.50 1
## [10,]  54.24  7.76 1
## [11,]  52.23  6.05 1
## [12,]  35.61  6.31 1
## [13,]  36.47  5.20 1
## [14,]  59.70 22.36 1
## [15,]  46.80 18.92 1
## [16,]  34.30  3.21 1
## [17,]  33.40  3.53 1
## [18,]  47.70  3.91 1
## [19,]  38.10 -3.65 1

#####################################
##   パッケージMASS ldaによる線形判別分析


suppressPackageStartupMessages(library(MASS))
dyall<-as.data.frame(yall)
dyall

##        x1    x2 z
## 1   15.50  1.14 0
## 2    5.60  1.03 0
## 3   28.30 -0.42 0
## 4   12.70 -8.05 0
## 5  -78.70  0.58 0
## 6   11.40  2.30 1
## 7   13.70  1.93 1
## 8   30.25  2.43 1
## 9   34.06  5.50 1
## 10  54.24  7.76 1
## 11  52.23  6.05 1
## 12  35.61  6.31 1
## 13  36.47  5.20 1
## 14  59.70 22.36 1
## 15  46.80 18.92 1
## 16  34.30  3.21 1
## 17  33.40  3.53 1
## 18  47.70  3.91 1
## 19  38.10 -3.65 1

(ZZZ<- lda(z~ . ,data=dyall))

## Call:
## lda(z ~ ., data = dyall)
##
## Prior probabilities of groups:
##       0       1
## 0.26316 0.73684
##
## Group means:
##       x1      x2
## 0 -3.320 -1.1440
## 1 37.711  6.1257
##
## Coefficients of linear discriminants:
##         LD1
## x1 0.033300
## x2 0.071224

(CC<-apply(ZZZ$means%*%ZZZ$scaling,2,mean))

##     LD1
## 0.75002#z=0.0333*(自己資本比率)+0.0712*(インタレストカバレ
#ッジ)-0.75

#ゼロが判別点となり、上式で得た得点がゼロ以下では倒産、
#ゼロより大きけば非倒産と判別される。定数項は倒産企業の
#各変数の平均値と非倒産企業の各説明変数の平均値との中点
#を式0.0333*(自己資本比率)+0.0712*(インタレストカバ
#レッジ)が通過するように計算されている。



##判別関数によるデータの判別結果
#左の1列目は当初の区分(倒産0 非倒産1) 2列目は
#モデルによる判別結果

data.frame(dyall[,3],predict(ZZZ)$class)

##    dyall...3. predict.ZZZ..class
## 1           0                  1
## 2           0                  1
## 3           0                  1
## 4           0                  0
## 5           0                  0
## 6           1                  1
## 7           1                  1
## 8           1                  1
## 9           1                  1
## 10          1                  1
## 11          1                  1
## 12          1                  1
## 13          1                  1
## 14          1                  1
## 15          1                  1
## 16          1                  1
## 17          1                  1
## 18          1                  1
## 19          1                  1

#度数分布表を見ると、判別点はゼロになるが若干この近辺で
# いくつかの企業が混在している
plot(ZZZ,dimen=1)

判別分析 プロビット分析  ロジット分析 ポアソン回帰

 

# 予測的中率をみるために下記のようなテーブルを作ってみ
#ると

table(実際=dyall[,3],予測=predict(ZZZ)$class)

##     予測
## 実際  0  1
##    0  2  3
##    1  0 14

#実際は倒産企業である会社をモデル予測では非倒産と判定し
#ている例が
3つあるが、非倒産企業についてはすべて的中さ
#せていることが分かる

###################
### Probit 分析の適用
##計算過程で警告が出るが計算手順の概略を示すことが目的な#のでそれを無視してすすめる。
suppressPackageStartupMessages(library(glm.predict))

## Error in library(glm.predict): there is no package #called 'glm.predict'

tousanprobit <- glm( z ~  x1 + x2,data = dyall, family = binomial(link = "probit"))


summary(tousanprobit)

 

##
## Call:
## glm(formula = z ~ x1 + x2, family = binomial(link = "probit"),
##     data = dyall)
##
## Deviance Residuals:
##     Min       1Q   Median       3Q      Max 
## -1.5934   0.0000   0.0026   0.0669   1.1268 
##
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -2.1174     1.8853   -1.12     0.26
## x1            0.1020     0.0824    1.24     0.22
## x2            0.4480     0.4313    1.04     0.30
##
## (Dispersion parameter for binomial family taken to be 1)
##
##     Null deviance: 21.9007  on 18  degrees of freedom
## Residual deviance:  7.8578  on 16  degrees of freedom
## AIC: 13.86
##
## Number of Fisher Scoring iterations: 11

#プロビットモデルによる各社の確率
#ここでは0.5を境界として0に近いほど倒産確率は高く、
#逆に1に近いほど非倒産の確率が高いと判断する


round(predict(tousanprobit,type="response"),digits=2)

##    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18   19
## 0.49 0.14 0.72 0.00 0.00 0.53 0.56 0.98 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.55

plot(fitted(tousanprobit))


table(実際 = dyall$z,
      予測 = round(fitted(tousanprobit)))

##     予測
## 実際  0  1
##    0  4  1
##    1  0 14

#この例では倒産企業の的中率がたまたま悪いが、実際の分析
#では
大量のデータを使うのでどの手法がよいかは分からな
#い。


#ロジットモデルの適用

tousanlogit <- glm(z  ~  x1 + x2,data = dyall, family = binomial(link="logit"))
summary(tousanlogit)

 

##
## Call:
## glm(formula = z ~ x1 + x2, family = binomial(link = "logit"),
##     data = dyall)
##
## Deviance Residuals:
##     Min       1Q   Median       3Q      Max 
## -1.6326  -0.0001   0.0284   0.1184   1.1052 
##
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -3.626      3.347   -1.08     0.28
## x1             0.176      0.147    1.20     0.23
## x2             0.779      0.745    1.05     0.30
##
## (Dispersion parameter for binomial family taken to be 1)
##
##     Null deviance: 21.9007  on 18  degrees of freedom
## Residual deviance:  7.9328  on 16  degrees of freedom
## AIC: 13.93
##
## Number of Fisher Scoring iterations: 9

 

#ここでも0.5を境界として0に近いほど倒産確率は高く、
#逆に1に近いほど非倒産の確率が高いと判断する


round(predict(tousanlogit,type="response"),digits=2)

##    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18   19
## 0.50 0.14 0.74 0.00 0.00 0.54 0.57 0.97 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.99 0.99 1.00 0.56

table(実際 = dyall$z,予測 = round(fitted(tousanlogit)))

##     予測
## 実際  0  1
##    0  4  1
##    1  0 14

#この例では倒産企業の的中率がよくなっているが、実際の分
#析では
大量のデータを使うのでなんともいえない。
#ロジットモデルの結果を3次元グラフで示すと以下のようにな
#る。ライブラリーlibrary(rockchalk)を使うと素人でもグラ
#フを描きやすい。


suppressPackageStartupMessages(library(rockchalk))

plotPlane(tousanlogit, plotx1 = "x1", plotx2 = "x2")

 

判別分析 プロビット分析  ロジット分析 ポアソン回帰

 

 

plotPlane(tousanlogit, plotx1 = "x1", plotx2 = "x2", phi = 20,theta=60)

 

判別分析 プロビット分析  ロジット分析 ポアソン回帰  

 

##ポアソン回帰の適用
#ポアソン回帰モデルは一定期間の事故発生件数、特許申請件
#数、犯罪発生件数など発生頻度の少ない事象をポアソン分布
#を仮定してモデル化されている。

#倒産についていえば一定期間に同一企業が何回も倒産するこ
#とはないのでカウントデータは0と1だけで2以上のケース
#は存在しないことになるが、幸にも,ポアソン分布の期待値
#(ここでは平均倒産回数)が小さい場合は回数が

#2以上の確率は無視できるくらい小さくなることが知られて
#おり、2値データ#モデルを十分近似できるといわれている。


tousanpo <- glm(z  ~  x1 + x2,data = dyall, family = poisson)
summary(tousanpo)

 

##
## Call:
## glm(formula = z ~ x1 + x2, family = poisson, data = dyall)
##
## Deviance Residuals:
##    Min      1Q  Median      3Q     Max 
## -1.143  -0.331  -0.013   0.274   0.718 
##
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.079844   0.702754   -1.54     0.12
## x1           0.023138   0.020934    1.11     0.27
## x2           0.000526   0.045798    0.01     0.99
##
## (Dispersion parameter for poisson family taken to be 1)
##
##     Null deviance: 8.5507  on 18  degrees of freedom
## Residual deviance: 5.7137  on 16  degrees of freedom
## AIC: 39.71
##
## Number of Fisher Scoring iterations: 5

round(fitted(tousanpo))

##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19
##  0  0  1  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1

table(実際 = dyall$z, 予測 = round(fitted(tousanpo)))

##     予測
## 実際  0  1
##    0  4  1
##    1  2 12

 

#ポアソン回帰で倒産か非倒産かの二者択一をさせている。

###########################################################
##ここからは全くの蛇足になるが、行列演算による判別分析の#係数等の推定値とライブラリによる係数のアウトプットの違
#いについて検討を加えてみた。

#多変量解析の教科書には線形代数を使った計算過程が示され
#ているので自己学習のために、上記と同じサンプルについて
#判別関数の係数を計算してみた。



#各郡の平均値の差 のベクトル
d1<-mean(ya1)-mean(yb1)
d1

## [1] -41.031

d2<-mean(ya2)-mean(yb2)
d2

## [1] -7.2697

D<-c(d1,d2)
D

## [1] -41.0314  -7.2697

# 平均偏差ベクトル
Sa1<-(ya1-mean(ya1))
Sa1

## [1]  18.82   8.92  31.62  16.02 -75.38

Sa2<-(ya2-mean(ya2))
Sa2

## [1]  2.284  2.174  0.724 -6.906  1.724

# 郡Aの行列
SSa<-cbind(Sa1,Sa2)
SSa

##         Sa1    Sa2
## [1,]  18.82  2.284
## [2,]   8.92  2.174
## [3,]  31.62  0.724
## [4,]  16.02 -6.906
## [5,] -75.38  1.724

SSa<-as.matrix(SSa,nrow(SSa),ncol(SSa))
SSa

##         Sa1    Sa2
## [1,]  18.82  2.284
## [2,]   8.92  2.174
## [3,]  31.62  0.724
## [4,]  16.02 -6.906
## [5,] -75.38  1.724

# 郡1の偏差平方和 行列
SSasqd<-t(SSa)%*%SSa
# 郡1の共分散行列  サンプル数-1 不偏推定量
acov<-SSasqd/(nrow(SSa)-1)
acov

##         Sa1     Sa2
## Sa1 1843.09 -38.830
## Sa2  -38.83  15.283

#郡2の偏差ベクトル
Sb1<-(yb1-mean(yb1))
Sb1

##  [1] -26.31143 -24.01143  -7.46143  -3.65143  16.52857  14.51857  -2.10143  -1.24143  21.98857
## [10]   9.08857  -3.41143  -4.31143   9.98857   0.38857

Sb2<-(yb2-mean(yb2))
Sb2

##  [1] -3.825714 -4.195714 -3.695714 -0.625714  1.634286 -0.075714  0.184286 -0.925714 16.234286
## [10] 12.794286 -2.915714 -2.595714 -2.215714 -9.775714

SSb<-cbind(Sb1,Sb2)
SSb

##             Sb1       Sb2
##  [1,] -26.31143 -3.825714
##  [2,] -24.01143 -4.195714
##  [3,]  -7.46143 -3.695714
##  [4,]  -3.65143 -0.625714
##  [5,]  16.52857  1.634286
##  [6,]  14.51857 -0.075714
##  [7,]  -2.10143  0.184286
##  [8,]  -1.24143 -0.925714
##  [9,]  21.98857 16.234286
## [10,]   9.08857 12.794286
## [11,]  -3.41143 -2.915714
## [12,]  -4.31143 -2.595714
## [13,]   9.98857 -2.215714
## [14,]   0.38857 -9.775714

# 郡2の行列
SSb<-as.matrix(SSb,nrow(SSb),ncol(SSb))
SSb

##             Sb1       Sb2
##  [1,] -26.31143 -3.825714
##  [2,] -24.01143 -4.195714
##  [3,]  -7.46143 -3.695714
##  [4,]  -3.65143 -0.625714
##  [5,]  16.52857  1.634286
##  [6,]  14.51857 -0.075714
##  [7,]  -2.10143  0.184286
##  [8,]  -1.24143 -0.925714
##  [9,]  21.98857 16.234286
## [10,]   9.08857 12.794286
## [11,]  -3.41143 -2.915714
## [12,]  -4.31143 -2.595714
## [13,]   9.98857 -2.215714
## [14,]   0.38857 -9.775714

# 郡2の偏差平方和 行列
SSbsqd<-t(SSb)%*%SSb
SSbsqd

##        Sb1    Sb2
## Sb1 2524.0 726.40
## Sb2  726.4 592.82

# 郡2の共分散行列  サンプル数-1 不偏推定量
bcov<-SSbsqd/(nrow(SSb)-1)
bcov

##         Sb1    Sb2
## Sb1 194.156 55.877
## Sb2  55.877 45.601

#郡1と郡2のプールされた分散行列 各郡の
#(サンプル数-1)での加重平均

Spl<-acov*(nrow(SSa)-1)/(nrow(SSa)+nrow(SSb)-2)+bcov*(nrow(SSb)-1)/(nrow(SSa)+nrow(SSb)-2)


Spl

##         Sa1    Sa2
## Sa1 582.141 33.593
## Sa2  33.593 38.468

#プールされた分散行列の逆行列
invSpl<-solve(Spl)


invSpl

##            Sa1        Sa2
## Sa1  0.0018090 -0.0015797
## Sa2 -0.0015797  0.0273755

#上記の逆行列に平均差ベクトルD を乗じ 係数を求める
(alpha<-invSpl%*%D)

##         [,1]
## Sa1 -0.06274
## Sa2 -0.13419#この係数はライブラリーMASSのldaを使って求めた数値
#x1 0.03329955 x2 0.07122402 と一見すると全く異なって
#いる。

# しかし、両者の係数ベクトルを作り正規化して比べると同一となっている。

# ldaで計算した係数ベクトルを作り正規化する


V1<-ZZZ$scaling
V1norm<-sqrt(t(V1)%*%V1)
(VV1<-V1/V1norm[1])

##        LD1
## x1 0.42353
## x2 0.90588

## 同様に行列計算で求めた係数を正規化すると
aa<-c(alpha[1:2])


(V<-aa)

## [1] -0.06274 -0.13419

Vnorm<-sqrt(t(V)%*%V)


(VV<-V/Vnorm[1])

## [1] -0.42353 -0.90588

##実質的に同じ結果を得る。正方行列の固有値λに対する固有
#ベクトル
のk倍もλに対する固有ベクトルとなるということ示
#している。

#念のためldaで求めた係数ベクトルを行列計算による係数ベク
#トルで除してみると、以下のように、すべて同一倍率になっ
#ている。


ZZZ$scaling/aa

##         LD1
## x1 -0.53076
## x2 -0.53076

       

金融工学の初等的計算例 目次

 

倒産確率 目次