/* 課題　10  解答例 */

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

/* (I) */
PDEex_put (S,K,r,sigma,T,M,N) :=
(
/* 2*N=原資産価格分割数 
     M=時間分割数　*/
/* 注: N>Mとする*/
Smax : 20*S                 /* 解くべき領域を原資産価格の 20 倍に設定 */,
xmax : log (Smax)        /*S_max に対応する x_max を計算*/,
dt : T/M                       /* 時間 [0,T] を M 個に分割 */,
dx : (xmax-log (S))/N    /* xmax から x=log (S) までを N 個に分割 */,
xmin : log (S) -N*dx      /* xmin か ら x=log (S) ま で を N 個 に 分 割 す る よ う に xmin を設定 */,
NN : 2*N,
x : makelist(xmin+i*dx, i, 0, NN)  /* 2N+1 次元ベクトル, x[1], x[2N+1] が境界，x[N+1]がlog(S) */,
dd : dt/ (dx^2)                            /* delta を定義する */,
a : (dd/2) * (sigma^2-dx* (r-sigma^2/2)),
b : (1-dd*sigma^2),
c : (dd/2) * (sigma^2+dx* (r-sigma^2/2)),

/* 初期条件 */
f : makelist(max (K-exp (x[i]),0), i, 1,NN+1)   /* t=T (満期)のペイオフを代入 */,
g : f                                                                  /*1 ステップ前の情報を記憶するために g を用意*/,

/*再帰計算*/
/* t->0 (j->M) に向かって解いていく*/
for j:1 thru M do (
f: makelist(a*g[i-1]+b*g[i]+c*g[i+1], i, 2, NN), 
g:f,
NN:NN-2     /* リストfの長さは，初項，最終項の２項について，２つ減少*/
),
P0 : exp (-r*T) *f[N-M+1] /* 原資産価格 S に対応するデリバティブ価格 */, 
/* BS 式計算 */
load(distrib),
d1 : (log (S/K) + (r+sigma^2/2) *T) / (sigma*sqrt (T)),
d2 : d1 - sigma*sqrt (T),
P0_A : exp (-r*T) *K*cdf_normal (-d2,0,1)-S*cdf_normal (-d1,0,1) ,
print("PDE=", P0,", ", "BS=", P0_A, ",", "Difference=", P0-P0_A)
);

/* Execution. */
PDEex_put(100,110,0.01,0.2,1,50,50);
PDEex_put(100,110,0.01,0.2,1,100,100);

kill(all);

/* (II) */
PDEex_American_put (S,K,r,sigma,T,M,N) :=
(
/* 2*N=原資産価格分割数 
     M=時間分割数　*/
/* 注: N>Mとする*/
Smax : 20*S                 /* 解くべき領域を原資産価格の 20 倍に設定 */,
xmax : log (Smax)        /*S_max に対応する x_max を計算*/,
dt : T/M                       /* 時間 [0,T] を M 個に分割 */,
dx : (xmax-log (S))/N    /* xmax から x=log (S) までを N 個に分割 */,
xmin : log (S) -N*dx      /* xmin か ら x=log (S) ま で を N 個 に 分 割 す る よ う に xmin を設定 */,
NN : 2*N,
x : makelist(xmin+i*dx, i, 0, NN)  /* 2N+1 次元ベクトル, x[1], x[2N+1] が境界，x[N+1]がlog(S) */,
dd : dt/ (dx^2)                            /* delta を定義する */,
a : (dd/2) * (sigma^2-dx* (r-sigma^2/2)),
b : (1-dd*sigma^2),
c : (dd/2) * (sigma^2+dx* (r-sigma^2/2)),
/* 初期条件 */
f : makelist( max (K-exp (x[i]),0), i, 1,NN+1)   /* t=T (満期)のペイオフを代入 */,
g : f                                                                  /*1 ステップ前の情報を記憶するために g を用意*/,
/*再帰計算*/
/* t->0 (j->M) に向かって解いていく*/
for j:1 thru M do (
f: makelist( max( exp(r*j*dt)*(K-exp(x[i])), a*g[i-1]+b*g[i]+c*g[i+1]), i, 2, NN),  
x: rest(rest(x,-1)), 
g:f,
NN:NN-2     /* リストfとxの長さは，初項，最終項の２項について，２つ減少*/
),
P0 : exp (-r*T) *f[N-M+1] /* 原資産価格 S に対応するデリバティブ価格 */, 
print("Price=", P0)
);

/* Execution. */
PDEex_American_put(100,110,0.01,0.2,1,50,50);
PDEex_American_put(100,110,0.01,0.2,1,100,100);

kill(all);

/* (III) */
PDEim_put(S,K,r,sigma,T,M,N):=
 (
Smax : 20*S
/* 解くべき領域を原資産価格の 20 倍に設定 */,
xmax : log (Smax)
/* Smax に対応する xmax を計算 */,
dt : T/M
/* 時間 [0,T] を M 個に分割 */,
dx : (xmax-log (S))/N
/* xmax から x=log (S) までを N 個に分割 */,
xmin : log (S) -N*dx
/* xmin か ら x=log (S) ま で を N 個 に 分 割 す る よ う に */
/* xmin を設定 */, 
NN : 2*N, 
x : makelist(xmin+i*dx,i, 0, NN), 

dd : dt/ (dx^2) /* delta を定義する */,
a : -dd* (sigma^2-dx* (r-sigma^2/2))/2,
b : (1+dd*sigma^2),
c : -dd* (sigma^2+dx* (r-sigma^2/2))/2,

/*  行列 D を生成 */
Da : addcol (addrow (zeromatrix(1, NN-2),diagmatrix(NN-2, a)),zeromatrix(NN-1,1)),
Db : diagmatrix (NN-1, b), 
Dc : addcol (zeromatrix(NN-1,1), addrow(diagmatrix(NN-2, c),zeromatrix(1, NN-2))),
D : Da + Db + Dc,
/* D の逆行列を計算 */
Dinv : invert(D), 
/* 初期条件 */
f0 : makelist(max (K- exp (x[i]),0), i, 1,NN+1) /* t=T (満期)のペイオフを代入 */,
f : rest(rest(f0,-1)) /*  境界 f0[1] (i=0),f0[NN+1] (i=2N) を除いた 2N-1 次元ベクトル */,
/*境界条件 */
g : makelist(0, NN-1) /* NN-1 個の要素0からなる数列を生成 */, 

/* t->0 (j->M) に向かって解いていく */
for j :1 thru M do (
g[1] : a*exp (r*j*dt) * f0[1] /* i=1 による境界条件を設定 */,
f : Dinv. (f-g)
),
P0 : exp (-r*T) *f[N][1],
/* 原資産価格 S に対応するデリバティブ価格 */
/* BS 式計算 */
load(distrib),
d1 : (log (S/K) + (r+sigma^2/2) *T) / (sigma*sqrt (T)),
d2 : d1 - sigma*sqrt (T),
P0_A : exp (-r*T) *K*cdf_normal (-d2,0,1)-S*cdf_normal (-d1,0,1) ,
print("PDE=", P0,", ", "BS=",P0_A, ",", "Difference=", P0-P0_A)
);

/* Execution */
PDEim_put (100,110,0.01,0.2,1,50,50);   
PDEim_put (100,110,0.01,0.2,1,100,100);

kill(all);

/* 課題IV */
PDEim_American_put(S,K,r,sigma,T,M,N) :=
(
Smax : 20*S  /*解くべき領域を原資産価格の20 倍に設定 */,
xmax : log (Smax)  /* S_max に対応するx_max を計算 */,
dt : T/M                /*時間[0,T] をM 個に分割 */,
dx : (xmax-log (S))/N   /*xmax からx=log (S) までをN 個に分割 */,
xmin : log (S) -N*dx    /*xmin からx=log (S) までをN 個に分割するようにxmin を設定 */,
NN  :2*N                     /* x=log (S) がメッシュ上に乗るために分割数をNN=2*N とする */,
x : makelist(xmin+i*dx ,i, 0, NN)  /*2N+1 次元ベクトル */,
dd : dt/ (dx^2)    /* delta  */, 
a : -dd* (sigma^2-dx* (r-sigma^2/2))/2,
b : (1+dd*sigma^2), 
c : -dd* (sigma^2+dx* (r-sigma^2/2))/2,

/*マトリクスD を生成*/
Da : addcol (addrow (zeromatrix(1, NN-2),diagmatrix(NN-2, a)),zeromatrix(NN-1,1)),
Db : diagmatrix (NN-1, b), 
Dc : addcol (zeromatrix(NN-1,1), addrow(diagmatrix(NN-2, c),zeromatrix(1, NN-2))),
D : Da + Db + Dc,

/* D の逆行列を計算 */
/* 関数: invert (M) あるは　M^^-1 行列 Mの逆行列を返す */
Dinv : invert(D), 
/* 初期条件 */
NNplus:NN+1, 
f0 : makelist(max (K- exp (x[i]),0), i, 1, NNplus) /*t=T の値，ペイオフを代入 */, 

f : rest(rest(f0,-1)) /*  境界 f0[1] (i=0),f0[NN+1] (i=2N) を除いた 2N-1 次元ベクトル */,
/*境界条件 */
NNminus:NN-1,
g : makelist(0, NNminus)  /* NN-1 個の要素0からなる数列を生成 */, 

/* t->0 (j->M) に向かって解いていく */
for j:1 thru M do(
/* 境界条件を設定*/
g[1] : a*exp (r*j*dt) * f0[1], 

for i:1 thru NNminus do (
f[i] : max (((Dinv . (f-g))[i])[1],exp(r*j*dt)*(K-exp(x[i]))) 
)), 
P0AP : exp (-r*T) *f[N] /* 原資産価格S に対応するデリバティブ価格 */, 
print ("a=", a,  "b=", b, "c=", c,"delta=", dd,"陰的有限差分=", P0AP)
);

/* Execution */
PDEim_American_put (100,110,0.01,0.2,1,50,50);
PDEim_American_put (100,110,0.01,0.2,1,100,100);