From 1c46773bb65074900da8de2e21cfbef5ca31c6bf Mon Sep 17 00:00:00 2001 From: Adrian Gao Date: Sat, 24 Jun 2023 18:06:43 +1000 Subject: [PATCH] WIP --- .../modified_merton/mod_merton_computation.py | 19 ++- .../mod_merton_create_lookup.py | 22 ++- .../modified_merton/mod_merton_simulation.py | 68 +++++++-- .../mod_single_cohort_computation.py | 133 ++++++++++++++++++ .../mod_single_cohort_solution.py | 15 ++ .../mod_single_merton_create_lookup.py | 49 +++++++ .../matlab/ModMertonComputation.m | 2 +- .../test_mod_merton_computation.py | 85 ++++++++++- 8 files changed, 358 insertions(+), 35 deletions(-) create mode 100644 src/frds/measures/modified_merton/mod_single_cohort_computation.py create mode 100644 src/frds/measures/modified_merton/mod_single_cohort_solution.py create mode 100644 src/frds/measures/modified_merton/mod_single_merton_create_lookup.py diff --git a/src/frds/measures/modified_merton/mod_merton_computation.py b/src/frds/measures/modified_merton/mod_merton_computation.py index 6af21d3b..928c2ef1 100644 --- a/src/frds/measures/modified_merton/mod_merton_computation.py +++ b/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,21 +229,21 @@ 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) @@ -251,7 +251,7 @@ def mod_merton_computation(fs, param, N, Nsim2, w=None, random_seed=1): 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() @@ -259,7 +259,6 @@ def mod_merton_computation(fs, param, N, Nsim2, w=None, random_seed=1): 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__": diff --git a/src/frds/measures/modified_merton/mod_merton_create_lookup.py b/src/frds/measures/modified_merton/mod_merton_create_lookup.py index 5d2707a0..c4a51440 100644 --- a/src/frds/measures/modified_merton/mod_merton_create_lookup.py +++ b/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 diff --git a/src/frds/measures/modified_merton/mod_merton_simulation.py b/src/frds/measures/modified_merton/mod_merton_simulation.py index 1a06b949..5e116889 100644 --- a/src/frds/measures/modified_merton/mod_merton_simulation.py +++ b/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__": diff --git a/src/frds/measures/modified_merton/mod_single_cohort_computation.py b/src/frds/measures/modified_merton/mod_single_cohort_computation.py new file mode 100644 index 00000000..2dd59680 --- /dev/null +++ b/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 diff --git a/src/frds/measures/modified_merton/mod_single_cohort_solution.py b/src/frds/measures/modified_merton/mod_single_cohort_solution.py new file mode 100644 index 00000000..e7c5ee34 --- /dev/null +++ b/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 diff --git a/src/frds/measures/modified_merton/mod_single_merton_create_lookup.py b/src/frds/measures/modified_merton/mod_single_merton_create_lookup.py new file mode 100644 index 00000000..efffaa05 --- /dev/null +++ b/src/frds/measures/modified_merton/mod_single_merton_create_lookup.py @@ -0,0 +1,49 @@ +import numpy as np +from numpy.random import RandomState + +from frds.measures.modified_merton.mod_single_cohort_computation import ( + mod_single_cohort_computation, +) + + +def mod_single_merton_create_lookup( + d, y, T, H, bookD, rho, ltv, xfs, xr, xF, xsig, N, Nsim2 +): + """ + same as `mod_merton_create_lookup` except that the function called differs + """ + # fmt: off + rng = RandomState(1) + w = rng.standard_normal((Nsim2, 3 * N)) + + J = xsig.shape[1] + K = xr.shape[2] + Q = xF.shape[3] + G = xfs.shape[0] + + fs = xfs[:, 0, 0, 0] + + xLt = np.zeros((G, J, K, Q)) + xBt = np.zeros((G, J, K, Q)) + xEt = np.zeros((G, J, K, Q)) + xFt = np.zeros((G, J, K, Q)) + xmdef = np.zeros((G, J, K, Q)) + xsigEt = np.zeros((G, J, K, Q)) + + 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] + + Lt, Bt, Et, _, _, _, sigEt, mFt, _, mdef, *_ = mod_single_cohort_computation(fs, param, N, Nsim2, w) + + 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 diff --git a/tests/measures/modified_merton/matlab/ModMertonComputation.m b/tests/measures/modified_merton/matlab/ModMertonComputation.m index 41d3c1f3..feccc9d4 100644 --- a/tests/measures/modified_merton/matlab/ModMertonComputation.m +++ b/tests/measures/modified_merton/matlab/ModMertonComputation.m @@ -1,5 +1,5 @@ -function [Lt, Bt, Et, LH, BH, EH, sigEt, mFt, def, mdef, face, FH, Gt, mu, F, sigLt] = ModMertonComputation(fs, param, N, Nsim2, w) +function [FHr2, Lt, Bt, Et, LH, BH, EH, sigEt, mFt, def, mdef, face, FH, Gt, mu, F, sigLt] = ModMertonComputation(fs, param, N, Nsim2, w) r = param(1); %log risk free rate diff --git a/tests/measures/modified_merton/test_mod_merton_computation.py b/tests/measures/modified_merton/test_mod_merton_computation.py index c95e3cec..b0e25dc2 100644 --- a/tests/measures/modified_merton/test_mod_merton_computation.py +++ b/tests/measures/modified_merton/test_mod_merton_computation.py @@ -1 +1,84 @@ -# This is the most difficult one... +import unittest +import pathlib +import numpy as np +from numpy.testing import assert_array_almost_equal + +try: + import matlab.engine + import matlab +except ImportError: + raise unittest.SkipTest("MATLAB Engine API not available.") + +from frds.measures.modified_merton.mod_merton_computation import mod_merton_computation + + +class ModMertonComputationTestCase(unittest.TestCase): + def setUp(self) -> None: + mp = pathlib.Path(__file__).parent.joinpath("matlab").as_posix() + self.eng = matlab.engine.start_matlab() + self.eng.cd(mp, nargout=0) + + def test_mod_merton_computation(self): + # fmt: off + fs = np.arange(-0.8, 0.85, 0.05) / (0.2 * np.sqrt(0.5) * np.sqrt(10)) + fs = fs.reshape(-1, 1) + + N = 10 # number of loan cohorts + Nsim2 = 10000 # number of simulated factor realization paths (10,000 works well) + + r = 0.01 # risk-free rate + d = 0.005 # depreciation rate of borrower assets + y = 0.002 # bank payout rate as a percent of face value of loan portfolio + T = 10 # loan maturity + bookD = 0.63 + H = 5 # bank debt maturity + D = bookD * np.exp(r * H) * np.mean(np.exp(r * np.arange(1, T+1))) # face value of bank debt: as if bookD issued at start of each cohort, so it grows just like RN asset value + rho = 0.5 + sig = 0.2 + ltv = 0.66 # initial ltv + g = 0.5 # prob of govt bailout (only used for bailout valuation analysis) + + j = 3 + bookF = 0.6+j*0.02 # initial loan amount = book value of assets (w/ coupon-bearing loans the initial amount would be book value) + + + param = [r, T, bookF, H, D, rho, ltv, sig, d, y, g] + + # Matlab results + result = self.eng.ModMertonComputation(fs, matlab.double(param), matlab.double(N), matlab.double(Nsim2), nargout=17) + r2, Lt_mt, Bt_mt, Et_mt, LH_mt, BH_mt, EH_mt, sigEt_mt, mFt_mt, default_mt, mdef_mt, face_mt, FH_mt, Gt_mt, mu_mt, F_mt, sigLt_mt = result + + # Python results + result = mod_merton_computation(fs, param, N, Nsim2) + r1, Lt, Bt, Et, LH, BH, EH, sigEt, mFt, default, mdef, face, FH, Gt, mu, F, sigLt = result + + n_precision = 9 + + # precision of 5 decimals + # assert_array_almost_equal(Lt, np.asarray(Lt_mt).ravel(), n_precision) + # assert_array_almost_equal(Bt, np.asarray(Bt_mt).ravel(), n_precision) + # assert_array_almost_equal(Et, np.asarray(Et_mt).ravel(), n_precision) + # assert_array_almost_equal(mFt, np.asarray(mFt_mt).ravel(), n_precision) + # assert_array_almost_equal(face, np.asarray(face_mt), n_precision) + # assert_array_almost_equal(FH, np.asarray(FH_mt), n_precision) + # assert_array_almost_equal(mu, np.asarray(mu_mt).ravel(), n_precision) # same + # assert_array_almost_equal(F, np.asarray(F_mt).ravel(), n_precision) + + assert_array_almost_equal(r1.flatten(), np.asarray(r2).flatten(), n_precision) # not equal to 5 decimals + print(r1.flatten()[:10]) + print(np.asarray(r2).flatten()[:10]) + # Hard to pass + # assert_array_almost_equal(LH, np.asarray(LH_mt), n_precision) # not equal to 5 decimals + # assert_array_almost_equal(BH, np.asarray(BH_mt), n_precision) # not equal to 5 decimals + # assert_array_almost_equal(EH, np.asarray(EH_mt), n_precision) # not equal to 5 decimals + # assert_array_almost_equal(sigEt, np.asarray(sigEt_mt).ravel(), n_precision) # not equal to 5 decimals + # assert_array_almost_equal(default, np.asarray(default_mt), n_precision) # FIXME: default is either 0 or 1, 0.00545% differ + # assert_array_almost_equal(mdef, np.asarray(mdef_mt).ravel(), n_precision) # not equal to 5 decimals + # assert_array_almost_equal(Gt, np.asarray(Gt_mt).ravel(), n_precision) # FIXME: very different + # assert_array_almost_equal(sigLt, np.asarray(sigLt_mt).ravel(), n_precision) # not equal to 5 decimals + + # fmt: on + + +if __name__ == "__main__": + unittest.main()