/* 課題7 */

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

/* (I)  */
mart_call2(r,u,S0,K,T,dt):=
(
d : 1/(1+u)-1,
q : (r*dt-d) / (u-d) /*リスク中立測度*/, 
N : truncate(T/dt) /* truncate () は引数の整数部分を返す関数 */, 
C0 :sum(max((1+u) ^j* (1+d) ^ (N-j) *S0-K,0) * binomial(N,j) *q^j* (1-q) ^ (N-j) / (1+r*dt)^N, j, 0, N),
print( " コールオプション価格=", C0 )
);

mart_call2(0.02,0.1,100,110,1,0.02);

/* (I) An Alternative */
mart_call_2a(r,u,S0,K,T,dt):=
(
d : 1/(1+u)-1,
N : truncate (T/dt), 
q : (r*dt -d) / (u-d), 
/* 危険資産の下方推移が行使価格を下回らない回数の最小値を求める*/
a : ceiling(log(K/S0/(1+d)^N )/ log((1+u)/(1+d))), 
/* そのときのコールオプションのペイオフを期待値へと加算していく*/
tmp : sum(
((1+u) ^j* (1+d) ^ (N-j) *S0-K) * binomial(N,j) *q^j* (1-q) ^ (N-j) / (1+r*dt)^N, j, a, N),
print(" マルチンゲール測度 q_u=",q, ", ", " コールオプション価格 =", tmp )
);

mart_call_2a(0.02,0.1,100,110,1,0.02);

/* (II)  */
Binomial_American_put(r,u,S,K,T,dt):=
(
d : 1/(1+u)-1,
N : truncate (T/dt), 
q : (r*dt -d) / (u-d), 
/* Calculation of the terminal values. */
value: makelist( max(K-(1+u)^ i * (1+d)^(N-i) * S , 0), i, 0, N ), 
/* backword recursive calculation of values. */
for n:1 thru N do 
 (
  m : (N-n) /* m denotes the time point from the initial time point 0.*/,
/* calculation of values at the time point m. */
  for i:1 thru (m+1) do 
   (
    value[i] : max( (q * value[i+1]+ (1-q) * value[i])/(1+r*dt), 
                            K-(1+u)^(i-1) * (1+d)^(m-i+1) * S)
   )
 ),
print("Price of the American put =",  value[1])
);

Binomial_American_put(0.02,0.1,100,110,1,0.02);

/* (II)  An Alternative */
Binomial_American_put_a(r,u,S,K,T,dt):=
(
d : 1/(1+u)-1,
N : truncate (T/dt), 
q : (r*dt -d) / (u-d), 
/* Calculation of the terminal values. */
value: makelist( max(K-(1+u)^ i * (1+d)^(N-i) * S , 0), i, 0, N ), 
/* backword recursive calculation of values. */
for n:1 thru N do 
 (
  m : (N-n) /* m denotes the time point from the initial time point 0.*/,
/* calculation of values at the time point m. */
  value: makelist(max( (q * value[i+1]+ (1-q) * value[i])/(1+r*dt), 
                            K-(1+u)^(i-1) * (1+d)^(m-i+1) * S), i, 1, m+1)
 ),
print("Price of the American put =",  value[1])
);

Binomial_American_put_a(0.02,0.1,100,110,1,0.02);

/* (III) */
Trinomial_American_put(r, sigma, S, K, T, dt):=
(
N : truncate(T/dt),
u : exp( sigma*sqrt( 3*dt ) ), 
d : 1/u,

/* calculation of risk neutral prob. */
qu : sqrt(dt/12)/sigma*(r-sigma^2/2)+1/6,
qn : 2/3,
qd : -sqrt(dt/12)/sigma*(r-sigma^2/2)+1/6,

/* 2*N+1 dim vector for the terminal values. */
value: makelist( max(K - u^(i-1) * d^N * S,0), i, 1, (2*N+1)),

/* backword recursive calculation of values. */
for n:1 thru N do
 (
  m : (N-n), 
/* m denotes the time point from the initial time point 0. */
/* calculation of values at the time point m.  */
 for i :1 thru (2*m+1) do 
   (
    value[i] : max((qu * value[i+2]
                   + qn * value[i+1]
                   + qd * value[i]
                   )/(1+r*dt), K - u^(i-1) * d^m * S)
   )
 ), 
print("Price of the American put price"=value[1])
);

Trinomial_American_put(0.02, 0.1, 100, 80, 1, 0.02);