WIP

src/frds/measures/modified_merton/mod_merton_computation.py

@@ -126,7 +126,7 @@ def mod_merton_computation(fs, param, N, Nsim2, w=None, random_seed=1):
 
     # From start of the first loan to the maturity of the first loan, first cohort only
     # Mingze's note: `xf1` unused, same in original Matlab code
-    xf1 = f[:, N] - f[:, 0]
+    # xf1 = f[:, N] - f[:, 0]
     # Factor shocks from the start of the first loan shocks until t
     f0w = np.tile(fw[:, N].reshape(fw[:, N].shape[0], 1), (1, N)) - fw[:, 0:N]
     # From the start of the first loan to the maturity of the first loan w/ shocks until t removed
@@ -229,37 +229,36 @@ def mod_merton_computation(fs, param, N, Nsim2, w=None, random_seed=1):
     # to one that depends only on (a)
     sc = L1 / bookF
     sc[ind1] = 1
-    FHr1 = FH1j
+    FHr1 = FH1j.copy()
     FHr1[ind2] = 0
     FHr2 = FH2j / sc
     FHr2[ind1] = 0
-    Lr1 = L1
+    Lr1 = L1.copy()
     Lr1[ind2] = 0
     Lr2 = L2 / sc
     Lr2[ind1] = 0
 
     LHj = np.zeros((Nsim2, N, Nsim1))
     for j in range(N):
-        FHr = np.reshape(FHr1[:, j, :] + FHr2[:, j, :], (Nsim2 * Nsim1, 1))
-        Lr = np.reshape(Lr1[:, j, :] + Lr2[:, j, :], (Nsim2 * Nsim1, 1))
+        FHr = np.reshape(FHr1[:, j, :] + FHr2[:, j, :], (Nsim2 * Nsim1, 1), order="F")
+        Lr = np.reshape(Lr1[:, j, :] + Lr2[:, j, :], (Nsim2 * Nsim1, 1), order="F")
         ind = np.argsort(FHr.flatten())
-        sortL = Lr[ind]
+        sortL = Lr[ind].copy()
         win = int(Nsim2 * Nsim1 / 20)  # / 10 seems to give about sufficient smoothness
         LHs = fftsmooth(sortL.flatten(), win)
         newInd = np.zeros(ind.shape, dtype=np.int64)
         for i in range(len(FHr)):
             newInd[ind[i]] = i
         LHsn = np.reshape(LHs[newInd], (Nsim2, Nsim1))
         LHsn = LHsn * sc[:, j, :]
-        LHj[:, j, :] = LHsn
+        LHj[:, j, :] = LHsn.copy()
 
     # Integrate over cohorts and discount to get portfolio payoff distribution at t+H
     LH1j = LHj.copy()
     LH1j[ind2] = 0
     LH2j = LHj.copy()
     LH2j[ind1] = 0
 
-    # FIXME: loss of precision here, LH1j is same as in Matlab
     LH = np.mean(
         LH1j * np.exp(-r * (rmat - HN) * (T / N)) + LH2j * np.exp(-r * (rmat - HN + N) * (T / N)),
         axis=1,
@@ -317,9 +316,9 @@ def mod_merton_computation(fs, param, N, Nsim2, w=None, random_seed=1):
     face = face[:, :szfs]
     Gt = Gt[:szfs]
 
-    # fmt: on
 
-    return Lt, Bt, Et, LH, BH, EH, sigEt, mFt, default, mdef, face, FH, Gt, mu, F, sigLt
+    return FHr2, Lt, Bt, Et, LH, BH, EH, sigEt, mFt, default, mdef, face, FH, Gt, mu, F, sigLt
+    # fmt: on
 
 
 if __name__ == "__main__":

src/frds/measures/modified_merton/mod_merton_create_lookup.py

@@ -14,7 +14,7 @@ def mod_merton_create_lookup(d, y, T, H, bookD, rho, ltv, xfs, xr, xF, xsig, N,
     Q = xF.shape[3]
     G = xfs.shape[0]
 
-    fs = xfs[:, 0, 0, 0] # (1,69,1,105), but in Matlab is (69,1,1,105)
+    fs = xfs[:, 0, 0, 0]
 
     xLt = np.zeros((G, J, K, Q))
     xBt = np.zeros((G, J, K, Q))
@@ -23,24 +23,20 @@ def mod_merton_create_lookup(d, y, T, H, bookD, rho, ltv, xfs, xr, xF, xsig, N,
     xmdef = np.zeros((G, J, K, Q))
     xsigEt = np.zeros((G, J, K, Q))
 
-    last_param = []
     for j in range(J):
         for k in range(K):
             for q in range(Q):
                 param = [xr[0, j, k, q], T, xF[0, j, k, q], H, bookD * np.exp(xr[0, j, k, q] * H), rho, ltv, xsig[0, j, k, q], d, y]
 
-                # FIXME: outputs dimenstions incorret: only 1*1, should be 69*1
                 Lt, Bt, Et, _, _, _, sigEt, mFt, _, mdef, *_ = mod_merton_computation(fs, param, N, Nsim2, w)
-                if param!=last_param:
-                    print(j,k,q, param, Et)
-                    last_param = param
-                xLt[:, j, k, q] = Lt.copy()
-                xBt[:, j, k, q] = Bt.copy()
-                xEt[:, j, k, q] = Et.copy() 
-                xFt[:, j, k, q] = mFt.copy()
-                xmdef[:, j, k, q] = mdef.copy()
-                xsigEt[:, j, k, q] = sigEt.copy()
 
-    return xLt, xBt, xEt, xFt, xmdef, xsigEt # FIXME: somehow all Lt, Bt, Et... remain same for all q
+                xLt[:, j, k, q] = Lt
+                xBt[:, j, k, q] = Bt
+                xEt[:, j, k, q] = Et
+                xFt[:, j, k, q] = mFt
+                xmdef[:, j, k, q] = mdef
+                xsigEt[:, j, k, q] = sigEt
+
+    return xLt, xBt, xEt, xFt, xmdef, xsigEt 
 
     # fmt: on

src/frds/measures/modified_merton/mod_merton_simulation.py

@@ -1,13 +1,18 @@
+import warnings  # fsolve may yield runtime warnings about invalid values
 import numpy as np
 from scipy.optimize import fsolve
 from scipy.interpolate import griddata
 from scipy.stats import norm
 
 from frds.measures.option_price import blsprice
+from frds.measures.modified_merton.merton_solution import merton_solution
 from frds.measures.modified_merton.mod_merton_computation import mod_merton_computation
 from frds.measures.modified_merton.mod_merton_create_lookup import (
     mod_merton_create_lookup,
 )
+from frds.measures.modified_merton.mod_single_merton_create_lookup import (
+    mod_single_merton_create_lookup,
+)
 
 
 def simulate():
@@ -74,18 +79,13 @@ def simulate():
     fs1 = []
     bookF1 = []
 
+    warnings.filterwarnings("ignore", category=RuntimeWarning)
+
     # Merton model fit to equity value and volatility from our model
     for k in range(K):
         initb = (Et[k] + D * np.exp(-r * H), sigEt[k] / 2)
 
-        def MertonSolution(b):
-            A, s = b
-            return [
-                (np.log(A) - np.log(D) + (r - y - s**2 / 2) * H) / (s * np.sqrt(H)),
-                norm.cdf(-mertdd[k]) - mertdef[k],
-            ]
-
-        bout = fsolve(MertonSolution, initb)
+        bout = fsolve(lambda b: merton_solution(b, Et[k], D, r, y, H, sigEt[k]), initb)
 
         A = bout[0]
         s = max(bout[1], 0.0001)
@@ -100,10 +100,14 @@ def MertonSolution(b):
         mertA[k] = A
         merts[k] = s
 
+    warnings.filterwarnings("default")
+
     mktCR = Et / (Et + Bt)
     sp = (1 / H) * np.log((D * np.exp(-r * H)) / Bt)
 
+    ###########################################################################
     # Alternative Model (1):
+    ###########################################################################
     # Perfectly correlated borrowers with overlapping cohorts
     rho = 0.99  # approx to rho = 1, code can't handle rho = 1
     sig = 0.20 * np.sqrt(0.5)
@@ -114,8 +118,9 @@ def MertonSolution(b):
     # so adjust bookD1 here
     bookD1 = bookD * np.mean(np.exp(r * np.arange(1, T + 1)))
 
-    # FIXME: problem is here, xfs,etc.dims don't match Matlab ndgrid's output
-    xfs, xsig, xr, xF = np.meshgrid(sfs, sig, r, np.arange(0.4, 1.45, 0.01))
+    xfs, xsig, xr, xF = np.meshgrid(
+        sfs, sig, r, np.arange(0.4, 1.45, 0.01), indexing="ij"
+    )
 
     xLt, xBt, xEt, xFt, xmdef, xsigEt = mod_merton_create_lookup(
         d, y, T, H, bookD1, rho, ltv, xfs, xr, xF, xsig, N, Nsim2
@@ -134,6 +139,49 @@ def MertonSolution(b):
 
     mdefsingle1 = ymdef
     mFtsingle1 = yFt
+    # recomputing Et, sigEt based on fs1, bookF1 gets back Et, sigEt to
+    # a very good approximation
+
+    ###########################################################################
+    # Model (2) Single cohort of borrowers
+    ###########################################################################
+    T = 5
+    rho = 0.5
+    sig = 0.2
+
+    xfs, xsig, xr, xF = np.meshgrid(
+        sfs, sig, r, np.arange(0.4, 1.45, 0.01), indexing="ij"
+    )
+
+    xLt, xBt, xEt, xFt, xmdef, xsigEt = mod_single_merton_create_lookup(
+        d, y, T, H, bookD1, rho, ltv, xfs, xr, xF, xsig, N, Nsim2
+    )
+
+    xxsigEt = xsigEt.squeeze()
+    xxEt = xEt.squeeze()
+    xxmdef = xmdef.squeeze()
+    xxFt = xFt.squeeze()
+
+    ymdef = griddata(
+        (xxEt.ravel(), xxsigEt.ravel()), xxmdef.ravel(), (Et, sigEt), method="cubic"
+    )
+    yFt = griddata(
+        (xxEt.ravel(), xxsigEt.ravel()), xxFt.ravel(), (Et, sigEt), method="cubic"
+    )
+
+    mdefsingle2 = ymdef
+    mFtsingle2 = yFt
+
+    ###########################################################################
+    # Model (3) Single (or perfectly correlated) borrower model
+    ###########################################################################
+
+    ###########################################################################
+    # Model (4) Merton model with asset value and asset volatility implied
+    # from modified model
+    ###########################################################################
+
+    pass
 
 
 if __name__ == "__main__":

src/frds/measures/modified_merton/mod_single_cohort_computation.py

@@ -0,0 +1,133 @@
+import numpy as np
+from scipy.optimize import fsolve
+
+from frds.measures.modified_merton.fftsmooth import fftsmooth
+from frds.measures.modified_merton.loan_payoff import loan_payoff
+from frds.measures.modified_merton.find_face_value_indiv import find_face_value_indiv
+
+
+def mod_single_cohort_computation(fs, param, N, Nsim2, w=None, random_seed=1):
+    # fmt: off
+    r = param[0]  # log risk-free rate
+    T = param[1]  # original maturity of bank loans
+    bookF = param[2]  # cash amount of loan issued = book value for a coupon-bearing loan issued at par
+    H = param[3]  # bank debt maturity
+    D = param[4]  # face value of bank debt
+    rho = param[5]  # borrower asset value correlation
+    ltv = param[6]  # initial LTV
+    sig = param[7]  # borrower asset value volatility
+    d = param[8]  # depreciation rate of borrower assets
+    y = param[9]  # bank payout rate
+
+    # With single borrower cohort, N just determines discretization horizon of
+    # shocks, not the number of cohorts 
+
+    # optional: calculate value of govt guarantee, with prob g, govt absorbs
+    # Loss given default (bank and equity values not adjusted for presence of
+    # govt guarantee) 
+    if len(param) > 10:
+        g = param[10]  # value of government guarantee
+    else:
+        g = 0
+
+    # optional: provide simulated factor shocks (for faster repeated computations) 
+    # if not provided as input, then generate here 
+    if w is None:
+        rng = np.random.RandomState(random_seed)
+        w = rng.normal(0, 1, (Nsim2, N))
+
+    # initial log asset value of borrower at origination
+    ival = np.log(bookF) - np.log(ltv)
+    sigf = np.sqrt(rho) * sig
+    # Mingze's note: `HN` must be an integer because it's used in array indexing
+    # This problem also presents in Nagel's original Matlab code
+    HN = H * (N / T) # maturity in N time
+    if isinstance(HN, float):
+        assert HN.is_integer()
+        HN = int(HN)
+    szfs = fs.shape[0]
+    fs = np.concatenate((fs, fs, fs), axis=0)
+    Nsim1 = fs.shape[0]
+
+    # Euler discretization 
+    f = np.concatenate(
+        (np.zeros((Nsim2, 1)), np.cumsum((r - d - 0.5 * sig ** 2) * (T / N) + sigf * np.sqrt(T / N) * w, axis=1)),
+        axis=1,
+    ) # FIXME:w is different 
+    # from start of first loan to maturity of first loan w/ shocks until t remove
+    f1 = f[:, N]
+
+    # add fs shocks after loan origination
+    fsa = (fs * sigf * np.sqrt(T)).reshape((1,Nsim1))
+    dstep = 10
+    df = sigf / dstep
+    fsa = fsa + df * np.concatenate(
+        (np.zeros((Nsim2, szfs)), np.ones((Nsim2, szfs)), -np.ones((Nsim2, szfs))),
+        axis=1,
+    )
+    f1j = np.tile(f1.reshape(Nsim2, 1), (1, Nsim1)) + fsa
+
+    initmu = r + 0.01
+    muout = fsolve(lambda mu: find_face_value_indiv(mu, bookF * np.exp(mu * T), ival, sig, T, r, d), initmu)
+    mu = muout[0]
+    F = bookF * np.exp(mu * T)
+
+    L1 = loan_payoff(F, f1j, ival, rho, sig, T)
+    face1 = np.ones_like(f1j) * F
+
+    ft = fsa + ival
+    fH1 = f[:, N+HN] 
+    fH1j = np.tile(fH1.reshape(Nsim2, 1), (1, Nsim1)) + fsa + ival
+    FH1j = np.exp(fH1j) * np.exp(0.5 * (1 - rho) * H * sig ** 2)
+    Ft = np.exp(ft)
+
+    FHr1 = FH1j
+    Lr1 = L1
+
+    LHj = np.zeros((Nsim2, Nsim1))
+    FHr = FHr1.flatten(order='F')
+    Lr = Lr1.flatten(order='F')
+    sort_indices = np.argsort(FHr)
+    sortF = FHr[sort_indices]
+    sortL = Lr[sort_indices]
+    win = (Nsim2 * Nsim1) // 20
+    LHs = fftsmooth(sortL, win)
+    new_indices = np.zeros_like(sort_indices)
+    new_indices[sort_indices] = np.arange(len(FHr))
+    LH1j = LHs[new_indices].reshape(Nsim2, Nsim1)
+
+    LH = LH1j * np.exp(-r * (T - H))
+    FH = FHr1
+    face = face1 * np.exp(-r * (T - H))
+
+    BH = np.minimum(D, LH * np.exp(-y * H))
+    EHex = LH * np.exp(-y * H) - BH
+    EH = LH - BH
+    GH = g * np.maximum(D - LH * np.exp(-y * H), 0)
+
+    Lt = np.mean(LH, axis=0) * np.exp(-r * H)
+    Bt = np.mean(BH, axis=0) * np.exp(-r * H)
+    Et = np.mean(EH, axis=0) * np.exp(-r * H)
+    Gt = np.mean(GH, axis=0) * np.exp(-r * H)
+    mFt = np.mean(Ft, axis=0)
+    def_val = np.ones_like(EHex)
+    ldef = EHex > 0
+    def_val[ldef] = 0
+    mdef = np.mean(def_val, axis=0)
+    sigEt = (dstep / 2) * (np.log(Et[szfs:2 * szfs]) - np.log(Et[2 * szfs:3 * szfs]))
+    sigLt = (dstep / 2) * (np.log(Lt[szfs:2 * szfs]) - np.log(Lt[2 * szfs:3 * szfs]))
+
+    Lt = Lt[:szfs]
+    Bt = Bt[:szfs]
+    Et = Et[:szfs]
+    LH = LH[:, :szfs]
+    BH = BH[:, :szfs]
+    EH = EH[:, :szfs]
+    FH = FH[:, :szfs]
+    mFt = mFt[:szfs]
+    def_val = def_val[:, :szfs]
+    mdef = mdef[:szfs]
+    face = face[:, :szfs]
+    Gt = Gt[:szfs]
+
+    return Lt, Bt, Et, LH, BH, EH, sigEt, mFt, def_val, mdef, face, FH, Gt, mu, F, sigLt

src/frds/measures/modified_merton/mod_single_cohort_solution.py

@@ -0,0 +1,15 @@
+from frds.measures.modified_merton.mod_single_cohort_computation import (
+    mod_single_cohort_computation,
+)
+
+
+def mod_single_cohort_solution(b, param, N, Nsim2, E, sigE):
+    fs = b[0]
+    param[2] = b[1]
+
+    # fmt: off
+    Lt, Bt, Et, LH, BH, EH, sigEt, mFt, default, mdef, face, FH, Gt, mu, F, sigLt = mod_single_cohort_computation(fs, param, N, Nsim2)
+    # fmt: on
+    err = [E - Et, sigE - sigEt]
+
+    return err