<jupyter_start><jupyter_text>Prj06**The question description**Approximate the price $e^{-rT}\mathbb E [F(S)]$, where- Asset follows $S = GBM(S_0, r, \sigma)$;- Knock-in payoff is given as$$F(S) = 1000 \cdot I(S_T > K) \cdot I\Big(\min_{1\le k \le m} S(t_k) < H\Big).$$__Parameters__- $r = 5\%, \sigma = 15\%, S(0) = 95$- $T = 0.25, m = 50, H = 85, K = 96.$- $k = 1000$, $n = 1000$__Todo__- Repeat $k$ times of the ordinary MC, with $n$ simulation for each MC. Find mean, MSE, and confidence interval using $k$ many MC outputs.- Repeat the above procedure with importance sampling.- Explain your decision of new probability in your importance sampling. Is it optimal choice? **Question1:**Repeat $k$ times of the ordinary MC, with $n$ simulation for each MC. Find mean, MSE, and confidence interval using $k$ many MC outputs. _____Analysis:__We denote the sequence of the payoff as ${X_1, X_2, \cdots, X_k}$. When we repeat k times of ordinary MC, by the CLT theorem, we know that\begin{equation*}\frac{\sqrt{k} (S_k - \mu)}{\sigma} \sim \mathcal{N}(0,1),\end{equation*}where $S_k = \frac{\sum_{i =1}^{k} X_i}{k}$ and $\mu = \mathbb{E}[X_i], \sigma = var[X_i].$ When the number of simulation is large enough, we can approximatively consider the distribution of the sample mean of payoff is normal with mean equal to the sample mean and variance equal to the sample varaince. So we get the bilateral confidence interval with confidence level equal to $1-\alpha$ as follow\begin{equation*}\mu = \frac{1}{k} \sum_{i =1}^{k} X_i \pm \frac{\sigma}{\sqrt{k}} Z(\alpha /2),\end{equation*}where $Z(\alpha /2)$ is the $\alpha/2$ upper quantile of standard normal distribution.<jupyter_code>import numpy as np
from scipy.stats import norm
import scipy.stats as ss
import matplotlib.pyplot as plt
from pylab import plt
plt.style.use('seaborn')
%matplotlib inline
## Exact simulation

# return: ndarray size (m, )
# Let W = bm_1d(T, m)
# W[i] is value at time (i+1)*dt
def bm_1d(T, m):
    dt = T/m
    dW = np.sqrt(dt)*norm.rvs(size = m)
    return dW.cumsum()

# arethmatic bm
def abm_1d(T, m, b, sigma, x0):
    dt = T/m
    dW = np.sqrt(dt)*norm.rvs(size = m)
    dX = b*dt + sigma*dW
    X = dX.cumsum()
    return x0 + X

#geometric bm
def gbm_1d(T, m, b, sigma, x0):
    L = abm_1d(T, m, b - 0.5*sigma**2, sigma, 0) #normalized log process
    return x0*np.exp(L)

## parameters
r = .05; sigma = 0.15; s0 = 95
T = 0.25; H = 85; K = 96
k = 1000; n = 1000; m = 50
dt = T/m
val = np.zeros(k)
payoff = np.zeros(n)
for i in np.arange(k):
    for j in np.arange(n):
        S = gbm_1d(T, m, r, sigma, s0)
        payoff[j] = 1000*float(S[-1]>K)*float(S.min()<H)
    val[i] = payoff.mean()*np.exp(-r*T)

print('MC computation for down_in_digital_call value: \n')
print('Mean is ' + str(val.mean()))    
print('MSE is ' + str(val.var()))

sample_mean = val.mean()
sample_variance = val.var()
## the plot of the value of knock-in call option in each simulation
plt.plot(np.linspace(0, k, k), val)
## the confidence level is 1-a, get the confidence interval of payoff
a = 0.05
Int_l = sample_mean - sample_variance/ np.sqrt(k) * ss.norm.cdf(a/2)
Int_r = sample_mean + sample_variance/ np.sqrt(k) * ss.norm.cdf(a/2)
print(f"The confidence interval of payoff with confidence level equal to 95% is:\n [{Int_l}, {Int_r}]")<jupyter_output>The confidence interval of payoff with confidence level equal to 95% is:
 [0.5609603571717193, 0.5787044245982197]<jupyter_text>**Question2:**Repeat the above procedure with importance sampling. ___**Analysis and mathematical dedection:**We denote $K = S_0 e^{c}$ and $H = S_0 e^{-b}$, the knock-in call option's payoff is given as$$\begin{align*} F(S) &=& 1000 \cdot I(S_T > S_0 e^{c}) \cdot I\Big(\min_{1\le k \le m} S(t_k) < S_0 e^{-b} \Big)  \\       &=& 1000 \cdot I(L_T > c) \cdot I\Big(\min_{1\le k \le m} L(t_k) < -b \Big) \\       &=& 1000 \cdot I(L_T > c) \cdot I(\tau \leq T),\end{align*} \tag{1}$$where $\tau = inf \{i : L(t_i) < -b \}.$We set that$$L_t = (r-\frac{1}{2} \sigma^2) t + \sigma W_t = \int_{0}^{t} (\hat{r_s} - 0.5 \sigma^2) ds + \sigma \hat{W_t},$$so we have $$W_t = \int_{0}^{t} \frac{\hat{r_s} - s}{\sigma} d s + \hat{W_t} =\int_{0}^{t} \theta_s d s + \hat{W_t}  .$$By the Girsanov theorem, we can get the Radon-Nikodym derivative as follow:$$\begin{equation*}\frac{d \mathbb{P}}{d \mathbb{\hat{P}}} = \gamma(\hat{W}) = exp \Big \{ \int_{0}^{T} \theta_s d \hat{W_s} - \frac{1}{2} \int_{0}^{T} \theta_{s}^{2} ds \Big \}.\end{equation*}$$And then we can get the price of the knock-in option by the formula$$e^{-rT} \hat{\mathbb{E}} \{ 1000 \cdot I(L_T > c) \cdot I(\tau \leq T) \gamma(\hat{W}) \}, $$where $L_t = \int_{0}^{t} (\hat{r_s} - 0.5 \sigma^2) ds + \sigma \hat{W_t}$ and $\tau = inf \{i : L(t_i) < -b \} \tag{2}$To get the exact form of the above formulas and to make the simulation easier to implement, we can take$$\hat{r_s} = \begin{cases}  -r_1, & s < h \\  r_2, & s \geq h\end{cases}$$then we have the form of $L_t$ as follow:$$L_t = \begin{cases}  (-r_1 - \frac{1}{2} \sigma^2) t + \sigma \hat{W_t}, & t < h \\  r_2(t -h) - r_1 h - \frac{1}{2} \sigma^2 t + \sigma \hat{W_t}, & t \geq h\end{cases}\tag{3}$$And since $$\theta_s = \begin{cases}  \theta_1 = \frac{-r_1 - r}{\sigma}, & s < h \\  \theta_2 = \frac{r_2 - r}{\sigma}, & s \geq h\end{cases}$$we can get the exact form of the $\gamma(\hat{W})$ by the formula$$\gamma(\hat{W}) = exp \Big \{ \theta_1 \hat{W_h} + \theta_2 (\hat{W_T} - \hat{W_h}) - \frac{1}{2} \theta_1^2 h - \frac{1}{2} \theta_2^2 (T-h) \Big \} \tag{4}$$ _____Algorithm:__Input $\leftarrow$ $m, n, k, S_0, r, \sigma, K, H, T, r_1, r_2, h$__Procedure__  Important sampling of MC  $\hspace{0.5cm}$ __for__ each $i \in (1,2,\cdots, k)$ __do__  $\hspace{1.0cm}$ __for__ each $j \in (1,2,\cdots, n)$ __do__  $\hspace{1.5cm}$ generate $\hat{W}$ on $(t_1 < t_2 < \cdots < t_m = T)$  $\hspace{1.5cm}$ generate $L$ on $(t_1 < t_2 < \cdots < t_m = T)$ by the formula $(3)$ $\hspace{1.5cm}$ get $\tau$ by the formula (2)  $\hspace{1.5cm}$ get the payoff by the formula (1)$\hspace{1.5cm}$ get $\gamma(\hat{W})$ by the formula (4)$\hspace{1.0cm}$ price estimator $\leftarrow$ $e^{-rT}$Average $(payoff \times \gamma)_{(1,2,\cdots, n)}$  $\hspace{0.5cm}$ mean $\leftarrow$ Average(price estimator) $\hspace{0.5cm}$ MSE $\leftarrow$ Variance(price estimator) __end Procedure__  Output mean, MSE___And we should realize that the value of $h, r_1, r_2$ may have huge influence to the result of the simulation. We should choose different value of $h, r_1, r_2$ to make the simulaiton more efficeint and reasonable.<jupyter_code>## get the parameter b and c in the formula

b = - np.log(H/s0)
c = np.log(K/s0)
print(f" The parameter b is equal to {b}")
print(f" The parameter c is equal to {c}")<jupyter_output>The parameter b is equal to 0.11122563511022437
 The parameter c is equal to 0.010471299867295437<jupyter_text>At first, we set $h = \frac{T}{2}$ set $r_1 = 0.06$ and $r_2 = 0.06$, and then we observe the result of the simulation<jupyter_code>## get the parameter theta_1 and theta_2 in the formula

r_1 = 0.06
r_2 = 0.06
theta_1 = - (r_1 + r)/sigma
theta_2 =  (r_2 - r) /sigma

h = int(m/2)<jupyter_output><empty_output><jupyter_text>The code of the question 2 is as follow:<jupyter_code>## importance sampling of MC simulation about the knoch-in call option with barrier H
## parameters
r = .05; sigma = 0.15; s0 = 95
T = 0.25; H = 85; K = 96
k = 1000; n = 1000; m = 50
dt = T/m
val = np.zeros(k)
payoff = np.zeros(n)
for i in np.arange(k):
    for j in np.arange(n):
        dW = np.sqrt(dt)*norm.rvs(size = m)
        W_t = dW.cumsum()      ## the path of Brownian Motion
        dL1_t = (- r_1 - 0.5 * sigma**2) *dt + sigma * dW[0:h]   ## the scatter of Lt in the first stage
        dL2_t = (r_2 - 0.5 *sigma**2)*dt + sigma* dW[h:(m+1)]       ## the scatter of Lt in the second stage
        L = np.append(dL1_t , dL2_t).cumsum()     ## the path of Lt
        gamma1 = theta_1 * dL1_t.cumsum()[-1] + theta_2 * dL2_t.cumsum()[-1] - (theta_1**2 + theta_2**2)/T/4
        gamma = np.exp(gamma1)
        payoff[j] = 1000*float(L[-1]>c)*float(L.min()<-b) * gamma 
    val[i] = payoff.mean()*np.exp(-r*T)
    
print('MC computation for down_in_digital_call value: \n')
print('Mean is ' + str(val.mean()))    
print('MSE is ' + str(val.var()))

sample_mean = val.mean()
sample_variance = val.var()
plt.plot(np.linspace(0, k, k), val)
## the confidence level is 1-a, get the confidence interval of payoff
a = 0.05
Int_l = sample_mean - sample_variance/ np.sqrt(k) * ss.norm.cdf(a/2)
Int_r = sample_mean + sample_variance/ np.sqrt(k) * ss.norm.cdf(a/2)
print(f"The confidence interval of payoff with confidence level equal to 95% is:\n [{Int_l}, {Int_r}]")<jupyter_output>The confidence interval of payoff with confidence level equal to 95% is:
 [0.5190260249280206, 0.5309204005470248]<jupyter_text>We can see that the MSE of the importance sampling of MC is 0.347 and it is less than the MSE of ordinary MC, so the importance sampling of MC can control the variate and more accuracy than the ordinary MC. And from the value plot of the knock-in option, we can see the scatter based on importance sampling of MC more centralized and less vibrate than the ordinary MC simulaiton. __Question3:__Explain your decision of new probability in your importance sampling. Is it optimal choice? ___As we related above, by the Girsanov theorem, we can get the Radon-Nikodym derivative as follow:$$\begin{align*}\frac{d \mathbb{P}}{d \mathbb{\hat{P}}} = \gamma(\hat{W}) = exp \Bigg \{ \int_{0}^{T} \theta_s d \hat{W_s} - \frac{1}{2} \int_{0}^{T} \theta_{s}^{2} ds \Bigg \}.\end{align*}$$And then we can get the price of the knock-in option by the formula$$e^{-rT} \hat{\mathbb{E}} \{ 1000 \cdot I(L_T > c) \cdot I(\tau \leq T) \gamma(\hat{W}) \}, $$where $L_t = \int_{0}^{t} (\hat{r_s} - 0.5 \sigma^2) ds + \sigma \hat{W_t}$ and $\tau = inf \{i : L(t_i) < -b \}.$And in the question 2 we set the $h = T/2$ and $$\hat{r_s} = \begin{cases}  -r_1, & s < h \\  r_2, & s \geq h\end{cases}$$which is a contant within the different two stage. So there may have some ameliorate to the simulation with a more desirable $\hat{r_s}.$ We can change the value of $r_1$ and $r_2$ to find the best parameters of the above function $\hat{r_s}$. More technically and reasonably, we can design $\hat{r_t}$ as a continous funciton and let it change with time.  What's more, we can change the time h to do some different simulation.<jupyter_code>## change the r_1 and r_2 and observe the change of MSE of the MC simulation 

r_list = np.linspace(0.03, 0.07, 5)

## importance sampling of MC simulation about the knoch-in call option with barrier H
## parameters
l = 5
val = np.zeros(k)
payoff = np.zeros(n)
sample_variance = np.zeros(l)
sample_mean = np.zeros(l)

for s in np.arange(l):
    for i in np.arange(k):
        for j in np.arange(n):
            dW = np.sqrt(dt)*norm.rvs(size = m)
            W_t = dW.cumsum()      ## the path of Brownian Motion
            dL1_t = (- r_list[s] - 0.5 * sigma**2) *dt + sigma * dW[0:h]   ## the scatter of Lt in the first stage
            dL2_t = (r_list[s] - 0.5 *sigma**2)*dt + sigma* dW[h:(m+1)]       ## the scatter of Lt in the second stage
            L = np.append(dL1_t , dL2_t).cumsum()     ## the path of Lt
            theta_1 = - (r_list[s] + r)/sigma
            theta_2 =  (r_list[s] - r) /sigma
            gamma1 = theta_1 * dL1_t.cumsum()[-1] + theta_2 * dL2_t.cumsum()[-1] - (theta_1**2 + theta_2**2)/T/4
            gamma = np.exp(gamma1)
            payoff[j] = 1000*float(L[-1]>c)*float(L.min()<-b) * gamma 
        val[i] = payoff.mean()*np.exp(-r*T)
    sample_mean[s] = val.mean()
    sample_variance[s] = val.var()

print(sample_mean)
print(sample_variance)<jupyter_output>[0.49518502 0.5119269  0.52221545 0.52986411 0.56765648]
[0.36998828 0.38018622 0.32946746 0.33276198 0.32219413]
{'fraction_non_alphanumeric': 0.169202319912583, 'fraction_numerical': 0.04757501891233084, 'mean_word_length': 3.4678933533608713, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 26, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '250706', 'license': 'MIT', 'n_tokens_mistral': 4949, 'n_tokens_neox': 4567, 'n_words': 1487, 'path': 'src/MA6628_JiaminJIAN_prj06.ipynb', 'repo_name': 'JiaminJIAN/20MA573', 'size': '11897'}