-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathq3.py
134 lines (117 loc) · 3.79 KB
/
q3.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Jan 1 15:18:30 2022
@author: benoitmuller
Question 3
"""
import numpy as np
import scipy as sp
import scipy.stats as st
import matplotlib.pyplot as plt
import time
import numpy.random as rnd
from mes_stats import RandomVariable
from sobol_new import *
from Payoff import Payoff
# File saving options:
save_figures = False # Change the value if needed
if (save_figures==True and
input("Do you really want to save figures"+
" into files?\n(yes/no): ") == "no"):
save_figures = False
alpha=0.01
NN=(2**np.arange(7,14)).astype(int) #(7,14)
mm=(2**np.arange(5,7)).astype(int) #(5,10)
MC1 = np.zeros((2,len(NN)))
MC2 = np.zeros((2,len(NN)))
QMC1 = np.zeros((2,len(NN)))
QMC2 = np.zeros((2,len(NN)))
plt.figure(figsize=(10, 9)) # figsize=(6,4) by default
plt.suptitle("Goal 3: Without pre-integration, estimated error" +
" with confidence $1-" + str(alpha) + "$")
for m in mm:
X=Payoff(m,K=500,S0=500)
for j in range(len(NN)):
MC1[:,j], MC2[:,j] = X.MC(NN[j],alpha)
QMC1[:,j], QMC2[:,j] = X.QMC(NN[j],alpha)
plt.subplot(221)
plt.loglog(NN,MC1[1,:],'.-',label='$m=$'+str(m))
plt.subplot(222)
plt.loglog(NN,QMC1[1,:],'.-',label='$m=$'+str(m))
plt.subplot(223)
plt.loglog(NN,MC2[1,:],'.-',label='$m=$'+str(m))
plt.subplot(224)
plt.loglog(NN,QMC2[1,:],'.-',label='$m=$'+str(m))
# Finest result (biggest m and N):
print("The computed intervals for m =", m,
", N =", NN[-1],"and alpha =",alpha,"are:")
print("V1: PIMC:",MC1[0,-1],"±",MC1[1,-1])
print(" PIQMC:",QMC1[0,-1],"±",QMC1[1,-1])
print("V2: PIMC:",MC2[0,-1],"±",MC2[1,-1])
print(" PIQMC:",QMC2[0,-1],"±",QMC2[1,-1])
# Error plots:
plt.subplot(221)
plt.loglog(NN,NN**(-0.5),"--",label='$1/\sqrt{N}$')
plt.title("Crude Monte Carlo for V1")
plt.xlabel('Sample size N')
plt.ylabel('Estimated error')
plt.legend()
plt.subplot(222)
plt.loglog(NN,NN**(-0.5),"--",label='$1/\sqrt{N}$')
plt.loglog(NN,10/NN,"--",label='$1/N$')
plt.title("Quasi Monte Carlo for V1")
plt.xlabel('Sample size N')
plt.ylabel('Estimated error')
plt.legend()
plt.subplot(223)
plt.loglog(NN,NN**(-0.5),"--",label='$1/\sqrt{N}$')
plt.title("Crude Monte Carlo for V2")
plt.xlabel('Sample size N')
plt.ylabel('Estimated error')
plt.legend()
plt.subplot(224)
plt.loglog(NN,NN**(-0.5),"--",label='$1/\sqrt{N}$')
plt.loglog(NN,10*1/NN,"--",label='$1/N$')
plt.title("Quasi Monte Carlo for V2")
plt.xlabel('Sample size N')
plt.ylabel('Estimated error')
plt.legend()
plt.tight_layout()
if save_figures == True:
plt.savefig('graphics/q3error.pdf')
# The interval plot:
plt.figure(figsize=(10,9))
plt.suptitle("Goal 1: Without pre-integration, interval of confidence $1-"
+ str(alpha) + "$, m="+str(mm[-1]))
plt.subplot(221)
plt.xscale('log')
plt.fill_between(NN,+MC1[0,:]+MC1[1,:],MC1[0,:]-MC1[1,:],alpha=0.5)
plt.plot(NN,MC1[0,:],'.-')
plt.title("Crude Monte Carlo for V1")
plt.xlabel('Sample size N')
plt.ylabel('mean and confidence interval')
plt.subplot(222)
plt.xscale('log')
plt.fill_between(NN,QMC1[0,:]+QMC1[1,:],QMC1[0,:]-QMC1[1,:],alpha=0.5)
plt.plot(NN,QMC1[0,:],'.-')
plt.title("Quasi Monte Carlo for V1")
plt.xlabel('Sample size N')
plt.ylabel('mean and confidence interval')
plt.subplot(223)
plt.xscale('log')
plt.fill_between(NN,+MC2[0,:]+MC2[1,:],MC2[0,:]-MC2[1,:],alpha=0.5)
plt.plot(NN,MC2[0,:],'.-')
plt.title("Crude Monte Carlo for V2")
plt.xlabel('Sample size N')
plt.ylabel('mean and confidence interval')
plt.subplot(224)
plt.xscale('log')
plt.fill_between(NN,QMC2[0,:]+QMC2[1,:],QMC2[0,:]-QMC2[1,:],alpha=0.5)
plt.plot(NN,QMC2[0,:],'.-')
plt.title("Quasi Monte Carlo for V2")
plt.xlabel('Sample size N')
plt.ylabel('mean and confidence interval')
plt.tight_layout()
if save_figures == True:
plt.savefig('graphics/q3interval.pdf')