# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python [default]
#     language: python
#     name: python3
# ---

# # HW 2: Corporate Bond Pricing (due by 9.21 Fri)

# We are going to compute the price of a corporate bond (subject to default) with Monte-Carlo simulation. Assume that 
# * the default time of a company follows the exponential distribution with intensity $\lambda=$__`def_rate`__. 
# * the riskfree interest rate is $r_f=$__`rf_rate`__ and the maturity of the bond is $T=$__`mat`__. 
# * in the case of default, you can recover some portion ($R=$__`recovery_rate`__) of the face value.
# * the coupon is 0%, i.e., it is a zero-coupon bond.
# * the face value of the bond is 1.0
# * use compound rate for discounting; the price of the default-free bond is $e^{-r_f T}$
#
# The Problem 1 of the [2017 ASP Midterm Exam](../files/ASP2017_Midterm.pdf) will be helpful.
#
# ### Instruction to upload your HW
# * Create a repository named __`PHBS_ASP_2018`__ (and clone it to your PC)
# * Copy this file to __`PHBS_ASP_2018/HW2/HW2.ipynb`__  (Please use the same name for repository and ipynb file)
# * Adding more code.
# * Run your your code to make sure that there's no error.
# * Upload (commit and sync) your file.

# ### 1. First, let's create a pricing function and check the std 

import numpy as np

def_rate = 0.1
rf_rate = 0.03
recovery = 0.3
mat = 10

#

# +
# First generate exponential random numbers
# Although you can generate directly using fault_time = np.random.exponential(scale=), let's use uniform random numbers.
n_sample = 10000
U = np.random.uniform(size=n_sample)
default_time = -(1/def_rate)*np.log(U)

# You can check if the RNs are correct by comparing the means
(default_time.mean(), 1/def_rate)


# +
# Put your code here to price the corporate bond

def corp_bond(def_rate=0.1, rf_rate=0.03, recovery=0.3, mat=10, n_sample=10000):
    U = np.random.uniform(size=n_sample)
    default_time = -(1/def_rate)*np.log(U)
    price = np.where(default_time>mat, 
                     np.exp(-rf_rate*mat), recovery*np.exp(-rf_rate*default_time))
    return np.mean(price)

# Call your function
corp_bond(def_rate, rf_rate, recovery, mat)

# Find the mean and std by calling the function 100 times.
arr = [corp_bond() for i in range(100)]
np.mean(arr), np.std(arr)
# -

# ### 2. Now, let's improve the function by reducing the MC variations.
# 1. Use antithetic method: If `U` is uniform random variable, so is `1-U`
# 2. Also shift the RNs to match the mean, `1/def_rate`

# +
# For example, antithetic method mean
n_sample = 10000
U = np.random.uniform(size=int(n_sample/2))
default_time = -(1/def_rate)*np.log(np.concatenate((U,1-U),axis=0))

# Mean-matching means
default_time += 1/def_rate-default_time.mean()
(default_time.mean(), 1/def_rate)


# +
# No include the two new features: `antithetic` and `mean_match`

def corp_bond_cv(def_rate=0.1, rf_rate=0.03, recovery=0.3, mat=10,
                 antithetic=False, mean_match=False, n_sample=10000):
    
    if(antithetic):
        U = np.random.uniform(size=int(n_sample/2))
        U = np.concatenate((U,1-U), axis=0)
    else:
        U = np.random.uniform(size=n_sample)

    default_time = -(1/def_rate)*np.log(U)
        
    if(mean_match):
        default_time += 1/def_rate - default_time.mean()
    
    price = np.where(default_time>mat, np.exp(-rf_rate*mat), recovery*np.exp(-rf_rate*default_time))
    return np.mean(price)        

# Find the mean and std by calling the function 100 times for (i) antithetic (ii) mean_match and (iii) both
arr1 = [corp_bond_cv(antithetic=True, mean_match=False) for i in range(100)]
arr2 = [corp_bond_cv(antithetic=False, mean_match=True) for i in range(100)]
arr3 = [corp_bond_cv(antithetic=True, mean_match=True) for i in range(100)]

np.array( [ np.mean(arr1), np.std(arr1), 
           np.mean(arr2), np.std(arr2), 
           np.mean(arr3), np.std(arr3) ] ).reshape([3,2])


# -

# ### 3. Finally, what is the analytic value of the corporate bond? How does it compare to your MC result above?

# +
### Put the analytic expression for the corporate bond price

def corp_bond_an(def_rate=0.1, rf_rate=0.03, recovery=0.3, mat=10):
    price = (recovery*def_rate*(1-np.exp((-rf_rate-def_rate)*mat)))/(rf_rate+def_rate)+np.exp((-rf_rate-def_rate)*mat)
    return np.mean(price)

corp_bond_an()
# -

#

#

StarcoderdataJupyterScriptsDedupFiltered