import numpy as np
import math
import cmUtilities as util
import numpy.linalg as anp
import importlib
import scipy
from scipy.stats import norm
import scipy.integrate as nInt
from scipy.stats import t as myT
import scipy.linalg as asp
from scipy.optimize import approx_fprime
import thresholdModels as th
import markovChain as mc
import time
import mixtureModels as mix
importlib.reload(util)
importlib.reload(th)
importlib.reload(mc)
importlib.reload(mix)
def jointDefaultProbability(p,q,myRho):
pr,err=nInt.quad(jointIntegrand,-5,5,args=(p,q,myRho))
return pr
def jointIntegrand(g,p,q,myRho):
p1 = th.computeP(p,myRho,g)
p2 = th.computeP(q,myRho,g)
f = p1*p2*util.gaussianDensity(g,0,1)
return f
def initializeRegion2r(N,rStart):
startRegion = np.zeros(N)
w = np.cumsum(rStart)
u = np.random.uniform(0,1,N)
for n in range(0,N):
if ((u[n]>0) & (u[n]<=w[0])):
startRegion[n] = 0
elif ((u[n]>w[0]) & (u[n]<1)):
startRegion[n] = 1
return startRegion
def createRatingData(K,N,T,Pin,wStart,myRho,nu,isT):
Q = cumulateTransitionMatrix(K,Pin)
if isT==1:
Delta = transformCumulativeTransitionMatrix_t(K,Q,nu)
else:
Delta = transformCumulativeTransitionMatrix(K,Q)
Y = np.zeros([N,T]) # latent variables
X = np.zeros([N,T]) # credit states
allP = np.zeros([N,T]) # default probabilities
Xlast = mc.initializeCounterparties(N,wStart) # initial states
X0 = Xlast
Plast = Pin[(Xlast-1).astype(int),-1]
for t in range(0,T):
Y[:,t] = th.getY(N,1,Plast,myRho,nu,isT)
for n in range(0,N):
if Xlast[n] == 4:
X[n,t] = 4
continue
else:
X[n,t] = migrateRating(Xlast[n],Delta,Y[n,t])
allP[:,t] = Pin[(Xlast-1).astype(int),-1]
Plast = allP[:,t]
Xlast = X[:,t]
return X,Y,Delta,allP,X0
def simulateCorrelatedTransitionData(K,N,T,Pin,wStart,myRho):
Q = cumulateTransitionMatrix(K,Pin)
Delta = transformCumulativeTransitionMatrix(K,Q)
Y = np.zeros([N,T]) # latent variables
X = np.zeros([N,T]) # credit states
allP = np.zeros([N,T]) # default probabilities
Xlast = mc.initializeCounterparties(N,wStart) # initial states
X0 = Xlast
Plast = Pin[(Xlast-1).astype(int),-1]
for t in range(0,T):
Y[:,t] = th.getY(N,1,Plast,myRho)
for n in range(0,N):
if Xlast[n] == 4:
X[n,t] = 4
continue
else:
X[n,t] = migrateRating(Xlast[n],Delta,Y[n,t])
allP[:,t] = Pin[(Xlast-1).astype(int),-1]
Plast = allP[:,t]
Xlast = X[:,t]
return X,Y,Delta,allP,X0
def createRatingData2r(K,N,T,P,wStart,rStart,myRho,nu,isT):
Q = cumulateTransitionMatrix(K,P)
Delta = transformCumulativeTransitionMatrix(K,Q)
rId = initializeRegion2r(N,rStart).astype(int)
Y = np.zeros([N,T]) # latent variables
X = np.zeros([N,T]) # credit states
allP = np.zeros([N,T]) # default probabilities
Xlast = mc.initializeCounterparties(N,wStart) # initial states
X0 = Xlast
Plast = P[(Xlast-1).astype(int),-1]
for t in range(0,T):
Y[:,t] = th.getY2r(N,1,Plast,myRho,rId,nu,P,isT)
for n in range(0,N):
if Xlast[n] == 4:
X[n,t] = 4
continue
else:
X[n,t] = migrateRating(Xlast[n],Delta,Y[n,t])
allP[:,t] = P[(Xlast-1).astype(int),-1]
Plast = allP[:,t]
Xlast = X[:,t]
return X,Y,Delta,allP,X0,rId
def migrateRating(lastX,Delta,myY):
transitionRow = (lastX-1).astype(int)
myMap = Delta[transitionRow,:]
if myY>=myMap[1]:
myX = 1
elif (myY<myMap[1]) & (myY>=myMap[2]):
myX = 2
elif (myY<myMap[2]) & (myY>=myMap[3]):
myX = 3
elif myY<myMap[3]:
myX = 4
return myX
def cumulateTransitionMatrix(K,M):
H = np.zeros([K,K])
for n in range(0,K):
for m in range(0,K):
H[m,(K-1)-n] = np.sum(M[m,(K-1)-n:K])
return H
def transformCumulativeTransitionMatrix(K,M_c):
H = np.zeros([K,K])
for n in range(0,K):
for m in range(0,K):
if M_c[n,m]>=0.9999999:
H[n,m]=5
elif M_c[n,m]<=0.0000001:
H[n,m] = -5
else:
H[n,m] = norm.ppf(M_c[n,m])
return H
def transformCumulativeTransitionMatrix_t(K,M_c,nu):
# Element-by-element inverse-normal transform
# of the cumulative transition matrix
H = np.zeros([K,K])
for n in range(0,K):
for m in range(0,K):
if M_c[n,m]>=0.9999999:
H[n,m] = 5
elif M_c[n,m]<=0.0000001:
H[n,m] = -5
else:
H[n,m] = myT.ppf(M_c[n,m],nu)
return H
def getSimpleEstimationData(T,X,allP):
N,T = X.shape
kVec = np.zeros(T)
nVec = np.zeros(T)
pVec = np.zeros(T)
kVec[0] = np.sum(X[:,0]==4)
nVec[0] = N
pVec[0] = np.mean(allP[:,0])
for t in range(1,T):
kVec[t] = np.sum(X[:,t]==4)-np.sum(X[:,t-1]==4)
nVec[t] = nVec[t-1] - kVec[t-1]
pVec[t] = np.mean(allP[X[:,t-1]!=4,t])
return pVec,nVec,kVec
def get2rEstimationData(T,X,X0,rId,allP,numP):
N,T = X.shape
kMat = np.zeros([T,numP])
nMat = np.zeros([T,numP])
pMat = np.zeros([T,numP])
for m in range(0,numP):
xLoc = (rId==m).astype(bool)
kMat[0,m] = np.sum(X[xLoc,0]==4)
nMat[0,m] = np.sum(xLoc)
pMat[0,m] = np.mean(allP[xLoc,0])
for t in range(1,T):
kMat[t,m] = np.sum(X[xLoc,t]==4)-np.sum(X[xLoc,t-1]==4)
nMat[t,m] = nMat[t-1,m] - kMat[t-1,m]
if np.sum(xLoc)==0:
pMat[t,m] = 0.0
else:
pMat[t,m] = np.mean(allP[(X[xLoc,t-1]!=4).astype(int),t])
return pMat,nMat,kMat
def getCMF(g,myRho,myP,myN,myK):
pg = th.computeP(myP,myRho,g)
f=util.getBC(myN,myK)*np.power(pg,myK)*np.power(1-pg,myN-myK)
cmf = f*util.gaussianDensity(g,0,1)
return cmf
def logLSimple(x,T,pVec,nVec,kVec):
L = np.zeros(T)
for t in range(0,T):
L[t],err = nInt.quad(getCMF,-5,5,
args=(x,pVec[t],nVec[t],kVec[t]))
logL = np.sum(np.log(L))
return -logL
def maxSimpleLogL(T,pVec,nVec,kVec):
myBounds = ((0.001,0.999),)
xStart = 0.5
r = scipy.optimize.minimize(logLSimple,
xStart,args=(T,pVec,nVec,kVec),
method='TNC',jac=None,bounds=myBounds,
options={'maxiter':1000})
return r.x,r.success
def computeSimpleScore(x0,T,pVec,nVec,kVec):
h = 0.00001
fUp = logLSimple(x0+h/2,T,pVec,nVec,kVec)
fDown = logLSimple(x0-h/2,T,pVec,nVec,kVec)
score = np.divide(fUp-fDown,h)
return score
def simpleFisherInformation(x0,T,pVec,nVec,kVec):
h = 0.000001
f = logLSimple(x0,T,pVec,nVec,kVec)
fUp = logLSimple(x0+h,T,pVec,nVec,kVec)
fDown = logLSimple(x0-h,T,pVec,nVec,kVec)
I = -np.divide(fUp-2*f+fDown,h**2)
return I
def getProdCMF(g,myRho,myP,myN,myK):
pg = th.computeP(myP,myRho,g)
return np.multiply(util.getBC(myN,myK),np.power(pg,myK)*np.power(1-pg,myN-myK))
def getCMF2r(g,myRho,pVec3,nVec3,kVec3):
myF=getProdCMF(g,myRho,pVec3,nVec3,kVec3)
return np.prod(myF)*util.gaussianDensity(g,0,1)
def logL2r(x,T,pMat,nMat,kMat):
L = np.zeros(T)
for t in range(0,T):
L[t],err = nInt.quad(getCMF2r,-5,5,args=(x,pMat[t,:],nMat[t,:],kMat[t,:]))
return -np.sum(np.log(L))
def log2rGridSearch(gridSize,numVar,T,pMat,nMat,kMat):
rhoRange = np.linspace(0.001,0.999,gridSize)
bigF = np.zeros([gridSize,gridSize])
rhoGrid = np.zeros([gridSize**numVar,numVar])
startTime = time.time()
for n in range(0,gridSize):
for m in range(0,gridSize):
print("Running: n:%d, m: %d" % (n+1,m+1))
rhoGrid = np.array([rhoRange[n],rhoRange[m]])
bigF[n,m] = -logL2r(rhoGrid,T,pMat,nMat,kMat)
print("Loop takes %d minutes." % ((time.time()-startTime)/60))
# Find the values associated with the biggest value of OF
bigMax = np.max(bigF)
for n in range(0,gridSize):
for m in range(0,gridSize):
if (bigF[n,m]==bigMax):
myN = n
myM = m
rhoStart = np.array([rhoRange[myN],rhoRange[myM]])
return rhoStart,bigF,rhoRange
def max2rLogL(T,pMat,nMat,kMat,xStart):
myBounds = ((0.001,0.999),(0.001,0.999))
r = scipy.optimize.minimize(logL2r,
xStart,args=(T,pMat,nMat,kMat),
method='TNC',jac=None,bounds=myBounds,
options={'maxiter':100})
return r.x,r.success
def hessian2r(x0,epsilon,T,pMat,nMat,kMat):
# The first derivative
f1 = approx_fprime(x0,logL2r,epsilon,T,pMat,nMat,kMat)
n = x0.shape[0]
hessian = np.zeros([n,n])
xx = x0
for j in range(0,n):
xx0 = xx[j] # Store old value
xx[j] = xx0 + epsilon # Perturb with finite difference
# Recalculate the partial derivatives for this new point
f2 = approx_fprime(xx, logL2r,epsilon,T,pMat,nMat,kMat)
hessian[:, j] = (f2 - f1)/epsilon # scale...
xx[j] = xx0 # Restore initial value of x0
return hessian
def score2r(x0,epsilon,T,pMat,nMat,kMat):
score = approx_fprime(x0,logL2r,epsilon,T,pMat,nMat,kMat)
return score
def mapRating(D,from_value,to_value,K):
if (to_value==K) & (from_value!=K):
d_u = D[from_value-1,to_value-1]
d_l = -5
elif (to_value==K) & (from_value==K):
d_u = -5
d_l = -5
else:
d_u = D[from_value-1,to_value]
d_l = D[from_value-1,to_value-1]
return d_l, d_u
def mapRatingData(Y,D,K):
N,T = Y.shape
d_low = np.zeros([N,T-1])
d_upp = np.zeros([N,T-1])
for n in range(0,N):
for m in range(1,T):
d_low[n,m-1],d_upp[n,m-1] = mapRating(D,
Y[n,m-1].astype(int),
Y[n,m].astype(int),K)
return d_low,d_upp
def buildCorrelationMatrix1F(x,N):
R = x*np.ones([N,N])+np.eye(N)
R[R==1+x] = 1
return R
def logLCopula(x,X,nu):
M,T = X.shape
R = buildCorrelationMatrix1F(x,M)
detR = anp.det(R)
Rinv = anp.inv(R)
V = 0
for t in range(0,T):
V += np.log(1+np.divide(np.dot(np.dot(X[:,t],Rinv),X[:,t]),nu))
return -(-0.5*(T*np.log(detR)+(nu+M)*V))
def scoreCopula(x0,myZ,nu):
h = 0.000000001
fUp = -logLCopula(x0+h/2,myZ,nu)
fDown = -logLCopula(x0-h/2,myZ,nu)
score = np.divide(fUp-fDown,h)
return score
def hessianCopula(x0,myZ,nu):
h = 0.0001
f = logLCopula(x0,myZ,nu)
fUp = logLCopula(x0+h,myZ,nu)
fDown = logLCopula(x0-h,myZ,nu)
I = np.divide(fUp-2*f+fDown,h**2)
return I
def maxCopulaLogL(startX,myZ,nu):
myBounds = ((0.001,0.999),)
r = scipy.optimize.minimize(logLCopula,
startX,args=(myZ,nu),
method='TNC',jac=None,bounds=myBounds,
options={'maxiter':1000})
return r.x,r.success
def mixtureMethodOfMoment(x,myP,myV,myModel):
if myModel==0: # Beta-binomial
M1 = mix.betaMoment(x[0],x[1],1)
M2 = mix.betaMoment(x[0],x[1],2)
elif myModel==1: # Logit
M1,err = nInt.quad(mix.logitProbitMoment,-8,8,args=(x[0],x[1],1,1))
M2,err = nInt.quad(mix.logitProbitMoment,-8,8,args=(x[0],x[1],2,1))
elif myModel==2: # Probit
M1,err = nInt.quad(mix.logitProbitMoment,-8,8,args=(x[0],x[1],1,0))
M2,err = nInt.quad(mix.logitProbitMoment,-8,8,args=(x[0],x[1],2,0))
elif myModel==3: # Poisson-gamma
M1 = mix.poissonGammaMoment(x[0],x[1],1)
M2 = mix.poissonGammaMoment(x[0],x[1],2)
elif myModel==4: # Poisson-lognormal
M1,err = nInt.quad(mix.poissonMixtureMoment,0.0001,0.9999,args=(x[0],x[1],1,0))
M2,err = nInt.quad(mix.poissonMixtureMoment,0.0001,0.9999,args=(x[0],x[1],2,0))
elif myModel==5: # Poisson-Weibull
M1,err = nInt.quad(mix.poissonMixtureMoment,0.0001,0.9999,args=(x[0],x[1],1,1))
M2,err = nInt.quad(mix.poissonMixtureMoment,0.0001,0.9999,args=(x[0],x[1],2,1))
f1 = M1 - myP
f2 =(M2-M1**2) - myV
return [1e4*f1, 1e4*f2]
def integrateGaussianMoment(g,r,myP,myMoment):
integrand = np.power(th.computeP(myP,r,g),myMoment)
return integrand*util.gaussianDensity(g,0,1)
def methodOfMomentsG(x,myP,myV):
if (x<=0) | (x>1):
return [100,100]
M1,err = nInt.quad(integrateGaussianMoment,-5,5,args=(x[0],myP,1))
M2,err = nInt.quad(integrateGaussianMoment,-5,5,args=(x[0],myP,2))
f1 = (M2-M1**2) - myV
return 1e4*f1
def thresholdMoment(g,w,p1,p2,myP,whichModel,myMoment,invCdf=0):
d1 = util.gaussianDensity(g,0,1)
if whichModel==1: # t
d2 = util.chi2Density(w,p2)
integrand = np.power(th.computeP_t(myP,p1,g,w,p2),myMoment)
if whichModel==2: # Variance-gamma
d2 = util.gammaDensity(w,p2,p2)
integrand = np.power(th.computeP_NVM(myP,p1,g,w,p2,invCdf),myMoment)
if whichModel==3: # Generalized hyperbolic
d2 = util.gigDensity(w,p2)
integrand = np.power(th.computeP_NVM(myP,p1,g,w,p2,invCdf),myMoment)
return integrand*d1*d2
def getThresholdMoments(x,myP,whichModel):
if whichModel==0: # Gaussian
M1,err = nInt.quad(integrateGaussianMoment,-5,5,args=(x[0],myP,1))
M2,err = nInt.quad(integrateGaussianMoment,-5,5,args=(x[0],myP,2))
elif whichModel==1: # t
lowerBound = np.maximum(x[1]-40,2)
support = [[-7,7],[lowerBound,x[1]+40]]
M1,err=nInt.nquad(thresholdMoment,support,args=(x[0],x[1],myP,whichModel,1))
M2,err=nInt.nquad(thresholdMoment,support,args=(x[0],x[1],myP,whichModel,2))
elif whichModel==2: # Variance-gamma
invCdf = th.nvmPpf(myP,x[1],0)
support = [[-7,7],[0,100]]
M1,err=nInt.nquad(thresholdMoment,support,args=(x[0],x[1],myP,whichModel,1,invCdf))
M2,err=nInt.nquad(thresholdMoment,support,args=(x[0],x[1],myP,whichModel,2,invCdf))
elif whichModel==3: # Generalized hyperbolic
invCdf = th.nvmPpf(myP,x[1],1)
support = [[-7,7],[0,100]]
M1,err=nInt.nquad(thresholdMoment,support,args=(x[0],x[1],myP,whichModel,1,invCdf))
M2,err=nInt.nquad(thresholdMoment,support,args=(x[0],x[1],myP,whichModel,2,invCdf))
return M1,M2
def thresholdMethodOfMoment(x,myP,myV,whichModel):
if (x[0]<=0) | (x[0]>1):
return 100
M1,M2 = getThresholdMoments(x,myP,whichModel)
f1 = M1 - myP
f2 =(M2-M1**2) - myV
return [1e4*f1,1e4*f2]
def getThresholdDefaultCorrelation(x,myP,whichModel):
if whichModel==0: # Gaussian
jp = th.jointDefaultProbability(myP,myP,x[0])
elif whichModel==1: # t
jp = th.jointDefaultProbabilityT(myP,myP,x[0],x[1])
elif whichModel==2: # Variance-gamma
jp = th.jointDefaultProbabilityNVM(myP,myP,x[0],x[1],0)
elif whichModel==3: # Generalized hyperbolic
jp = th.jointDefaultProbabilityNVM(myP,myP,x[0],x[1],1)
return np.divide(jp-myP**2,myP*(1-myP))
def getCMF_cr(s,myW,myA,myP,myN,myK):
ps=myP*(1-myW)+myP*myW*s
f=util.getBC(myN,myK)*np.power(ps,myK)*np.power(1-ps,myN-myK)
return f*util.gammaDensity(s,myA,myA)
def logLSimple_cr(x,T,pVec,nVec,kVec,myA):
L = np.zeros(T)
for t in range(0,T):
L[t],err = nInt.quad(getCMF_cr,0,100,
args=(x,myA,pVec[t],nVec[t],kVec[t]))
return -np.sum(np.log(L))
def crPlusMoment(s,myW,myA,myP,momentNumber):
v0 = myP*(1-myW) + myP*myW*s
myDensity = util.gammaDensity(s,myA,myA)
return np.power(v0,momentNumber)*myDensity
def crPlusScore(x0,T,pVec,nVec,kVec,myA):
h = 0.0000001
fUp = logLSimple_cr(x0+h/2,T,pVec,nVec,kVec,myA)
fDown = logLSimple_cr(x0-h/2,T,pVec,nVec,kVec,myA)
return np.divide(fUp-fDown,h)
def crPlusFisherInformation(x0,T,pVec,nVec,kVec,myA):
h = 0.000001
f = logLSimple_cr(x0,T,pVec,nVec,kVec,myA)
fUp = logLSimple_cr(x0+h,T,pVec,nVec,kVec,myA)
fDown = logLSimple_cr(x0-h,T,pVec,nVec,kVec,myA)
return -np.divide(fUp-2*f+fDown,h**2)
def jointDefaultProbabilityRegion(p,q,rhoVec):
pr,err=nInt.quad(jointIntegrandRegion,-10,10,args=(p,q,rhoVec))
return pr
def jointIntegrandRegion(g,p,q,rhoVec):
p1 = th.computeP(p,rhoVec[0],g)
p2 = th.computeP(q,rhoVec[1],g)
density = util.gaussianDensity(g,0,1)
f = p1*p2*density
return f
def getRegionalDefaultCorrelation(rhoVec,myP):
jp = jointDefaultProbabilityRegion(myP[0],myP[1],rhoVec)
return np.divide(jp-myP[0]*myP[1],np.sqrt(myP[0]*(1-myP[0]))*np.sqrt(myP[1]*(1-myP[1])))
