/* 課題　9  解答例 */

/* ［数値処理］→［自動的に数値で出力］*/
if numer#false then numer:false else numer:true;
/*浮動小数点の桁数=7に指定 */
fpprec : 7;
/* 浮動小数点表示桁数を7に指定*/
fpprintprec:7;

/* (I) */
put_monte_BS(S,K,r,sigma,T,n):=
(
/* generation of n normal random numbers */
load(distrib),
x : random_normal (0,1, n) , 
/* standarization */
y : (x - mean (x))/std (x) , 
/* caluculation of paypffs */
z : K - S*exp ((r-0.5*sigma^2) *T+sigma*sqrt (T) *y) , 
/* option price by the Monte-Carlo */
M_P0 : sum (if z[i]>0 then z[i] else 0, i, 1, n) *exp (-r*T) /n, 
/*  option price by the BS formula */
d1 : (log(S/K)+(r+sigma^2/2)*T) / (sigma*sqrt(T)),
d2 : d1 - sigma * sqrt(T),
BS_P0 : exp (-r*T) * K * cdf_normal(-d2, 0,1 ) - S * cdf_normal(-d1, 0, 1), 
print("Price by the Monte-Carlo=", M_P0, ", " , "Price by the BS formula=", BS_P0)
);

/* Execution */
put_monte_BS(100,100,0.01,0.2,1,10000);

kill(all);

/* (II) */
rainbow_put_option(S,K,r,Sigma,T,n):=
(
/* S = vector of the current underlyings prices
     Sigma = matrix of the underlyings volatilities
     m = number of the underlyings */
m  : length (S),
/* calculation of each underlyings volatility by the inner product */
v : makelist((Sigma.transpose(Sigma))[i, i], i, 1, m), 

payoff:0,
load(distrib),
for i:1 thru n do (
/*  generation of m normal random numbers */
epsilon : random_normal (0,1,m), 
/* maximum of m underlyings prices at the terminal */
tmp : lmax(S*exp ((r-0.5*v) *T+(epsilon. transpose(Sigma) )[1]*sqrt (T))),
payoff: payoff +max(K-tmp,0)
),
price : payoff/n*exp (-r*T), 
print ("price=", price)
);

/* Execution */
Sigma : matrix ([0.2,0.1],[0.4,0.5]);
rainbow_put_option ([100,100],100,0.01,Sigma,1,10000);

kill(all);

/* (III) */
basket_put_option(S,K,r,Sigma,T,n):=
(
/* S = vector of the current underlyings prices
     Sigma = matrix of the underlyings volatilities
     m = number of the underlyings */
m  : length (S),
/* calculation of each underlyings volatility by the inner product */
v : makelist((Sigma.transpose(Sigma))[i, i], i, 1, m), 

load(distrib),
payoff:0, 
for i:1 thru n do
(
 /* summation of the underlyings prices at the terminal */
 epsilon : random_normal(0,1, m),
 tmp : S . exp ((r-0.5*v) *T+(epsilon.transpose(Sigma))[1]*sqrt (T)),
/* sequential summation of the payoffs */
 payoff : payoff + max (K-tmp/m,0)
),
price : payoff/n*exp (-r*T),
print("price=", price)
);

/*  Execution */
Sigma : matrix ([0.2,0.1],[0.4,0.5]);
basket_put_option ([100,100],100,0.01,Sigma,1,10000);