信用リスクモデリング　－測定と管理ー　森平爽一郎先生

#p.32 2.2.2 実績デフォルト相関

data(anscombe)

Y=as.matrix(round(anscombe[,!(colnames(anscombe) %in% c("x1","x2","x3","x4"))]))

rownames(Y)=paste0("t",c(1:nrow(Y)))

N=apply(Y,1,max)+1;FD=Y/N

cor_mat=array(0,dim=c(ncol(Y),ncol(Y)))

for(i in 1:ncol(cor_mat)){
for(j in 1:nrow(cor_mat)){  

cor_mat[i,j]=(mean(FD[,i]*FD[,j])-mean(FD[,i])*mean(FD[,j]))/(sqrt(mean(FD[,i])*(1-mean(FD[,i])))*sqrt(mean(FD[,j])*(1-mean(FD[,j]))))  

}}

#信用リスクモデリング　－測定と管理ー　森平爽一郎著

#p.40 2.4.2 「非線形実績デフォルト率」推定モデル

#重み付き最小二乗法による反復解法

data(airquality)

airquality[is.na(airquality)>0]=0

x=round(airquality[,colnames(airquality) %in% c("Ozone")])

N=max(x)+10

D=x;FD=D/N+1/(2*N)

y=as.matrix(cbind(rep(1,nrow(airquality)),airquality[,!(colnames(airquality) %in% c("Ozone"))]))

odds=log(FD/(1-FD))

V_ep=rep(0,nrow(y))

ite=20

for(l in 1:ite){

V_ep_pre=V_ep

beta=solve(t(y)%*%y)%*%t(y)%*%odds

pre=y%*%beta

ep=pre-log(FD/(1-FD))

prob=1/(1+exp(-pre))

#taylor展開での近似

ep_approx=(prob-FD)/(FD*(1-FD))

V_ep=1/(N*prob*(1-prob))

#重みを計算

w=sqrt(V_ep)

odds=odds*w;y=y*c(w)

print(sum(abs(V_ep_pre-V_ep)))

}

#p.48 3.1.2

#倒産・非倒産

Y=c(rep(1,7),rep(0,7))

#負債比率

X=c(1.932,0.993,1.053,0.881,0.896,0.693,0.588,0.403,0.617,0.563,0.264,0.55,0.72,0.67)

a=cov(X,Y)/var(X)

b=mean(Y)-mean(X)*a

PD_hat=b+a*X

#誤差項の分散はベルヌーイ分布のため、異なる（不均一分散）

V_ep=PD_hat*(1-PD_hat)

#加重最小二乗法（均一分散にするため） 微妙な出来。

X=X[-1];Y=Y[-1]

PD_hat_pre=rep(0,length(X))

ite=10

for(l in 1:ite){

a=cov(X,Y)/var(X)

b=mean(Y)-mean(X)*a

PD_hat=b+a*X

sd_ep=sqrt(PD_hat*(1-PD_hat))

sd_ep[is.na(sd_ep)>0]=1

X=X/sd_ep

Y=Y/sd_ep

print(sum(abs(PD_hat-PD_hat_pre)))

PD_hat_pre=PD_hat 

}

#倒産・非倒産

Y=c(rep(1,7),rep(0,7))

#負債比率

X=c(1.932,0.993,1.053,0.881,0.896,0.693,0.588,0.403,0.617,0.563,0.264,0.55,0.72,0.67)

X=as.matrix(cbind(rep(1,length(Y)),X))

#bernoulli-logistic regression p.105

#one dimension

alpha=1;beta=1

ite=10^2

eta=0.1

for(j in 1:ite){

prob=1/(1+exp(-(X%*%c(alpha,beta))))

V=diag(c(prob*(1-prob)))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%(Y-prob)

alpha=alpha+eta*mat[1];beta=beta+eta*mat[2]

loglik=sum(Y*log(prob)+(1-Y)*log(1-prob))

print(loglik)

}

PD=prob

#p.57

#オッズ比に対する感応度

beta

#k番目のリスクファクター（特徴量）が微小変化した際の影響

delta=beta*PD*(1-PD)

#デルタの信用リスクファクターの変化に対する影響（ガンマ）

gamma=beta*PD*(1-PD)*(1-2*PD)

#6.4.4 最尤法による推定　p.128

data("airmiles")

Et=airmiles

len=length(airmiles)

ts=50

t=len

Dt=3000

n=15

times=c(1:n)

A=Et[1:n]+Dt

A_vec=c()

for(j in 1:(len-n)){

A_pre=A

while(sum(abs(A-A_pre))>10^(-3))({ 

  A_pre=A

  rA=A[2:length(A)]/A[1:(length(A)-1)]-1

  mu_A=mean(rA);sig_A=sd(rA)

  d1=(log(A/Dt)+(mu_A+sig_A^2/2)*(t-times))/(sig_A*sqrt(t-times))

  d2=d1-sig_A*sqrt(t-times)

  A=(Et[1:n]-Dt*exp(-mu_A*(t-times))*pnorm(d2))/pnorm(d1)

  PD=1-pnorm(d2)

  print(sum(abs(A-A_pre)))

})

Dt=A[1]/exp(sig_A*sqrt(ts-times[1])*d2[1]-(mu_A-sig_A^2/2)*(ts-times[1]))

A=c(A,Et[n+j]+Dt)

times=c(times,max(times)+1)

}

rA=A[2:length(A)]/A[1:(length(A)-1)]-1

mu_A=mean(rA);sig_A=sd(rA)

d1=(log(A/Dt)+(mu_A+sig_A^2/2)*(t-times))/(sig_A*sqrt(t-times))

d2=d1-sig_A*sqrt(t-times)

PD=1-pnorm(d2)