Hi, not wanting to start the new year with doubt, I have rewritten along the lines of the VI-adap.edp example the variational inequality of the basket price of American options. Given my inexperience, would anyone be able to tell me if this is a good idea?
The results obtained are similar to the code found on the internet which had a slightly different implementation from the theory. Could you please check to see if I have written any corruptions?
Thank you and happy new year.
//--- algo parameters
int m=40;
int L=80;
int LL=80;
int j=-1;
int kmax=1;
int k=0;
//--- financial parameters
real T=1;
real sigmax=0.35;
real sigmay=0.3;
real rho=-0.3;
real r=0.02;
real K=40;
real dt=0.01;
//---
bool debug = false;
//mesh Th=square(20,20);
mesh Th=square(m,m,[L*x,LL*y]);
real eps=0.35;
fespace Vh(Th,P1); // P1 FE space
int n = Vh.ndof; // number of Degree of freedom
Vh uh,uhp,vh; // solution and previous one
Vh Ik; // to def the set where the containt is reached.
real[int] rhs(n); // to store the right and side of the equation
real c=1000; // the parameter of the algoritme
func f=1; // right hand side function
func fd=0; // Dirichlet boundary condition function
Vh g=1;
// array to store
real[int] Aii(n),Aiin(n); // store the diagonal of the matrix
real tol=0.001,tolmin=0.0001;
real tgv = 1e30; // a huge value of exact penalisation of boundary condition
real res=0.;
// American option early exercise and other
func e = max(K-max(x,y),0.);
Vh pay=e;
uh=pay;
Vh xveloc = -x*r+x*sigmax^2+x*rho*sigmax*sigmay/2;
Vh yveloc = -y*r+y*sigmay^2+y*rho*sigmax*sigmay/2;
// the variational form of the problem:
varf a(uh,vh) = int2d(Th)( uh*vh*(r+1/dt)
+ dx(uh)*dx(vh)*(x*sigmax)^2/2 + dy(uh)*dy(vh)*(y*sigmay)^2/2
+ (dy(uh)*dx(vh) + dx(uh)*dy(vh))*rho*sigmax*sigmay*x*y/2)
- int2d(Th)(vh*convect([xveloc,yveloc],dt,uhp)/dt)
+ on(2,3,uh=0);
// two version of the problem
matrix A=a(Vh,Vh,tgv=tgv,solver=CG);
matrix AA=a(Vh,Vh);
// the mass Matrix construction:
varf vM(uh,vh) = int2d(Th)(uh*vh);
matrix M=vM(Vh,Vh); // to do a fast computing of $L^2$ norm : sqrt( u'*(w=M*u))
Aii=A.diag; // get the diagonal of the matrix
rhs = a(0,Vh,tgv=tgv);
Ik =0;
uhp=0;
Vh lambda=0;
int kadapt=0,kkadapt=0, iter=0;
real[int] b(n) ;
real[int] Ak(n); // the complementary of Ik ( !Ik = (Ik-1))
j=0;
// iterate for all time T
// Start counting time
real time=clock();
while (iter*dt <= T)
{
k=0;
while(k<kmax)
{
rhs = a(0,Vh,tgv=tgv);
b=rhs; // get a copy of the Right hand side
// The operator Ik- 1. is not implement so we do:
Ak= 1.; Ak -= Ik[]; // build Ak = ! Ik
// adding new locking condition on b and on the diagonal if (Ik ==1 )
b = Ik[] .*pay[];
b *= tgv;
b -= Ak .* rhs;
Aiin = Ik[] * tgv;
Aiin += Ak .* Aii; //set Aii= tgv $ i \in Ik $
A.diag = Aiin; // set the matrix diagonal
set(A,solver=CG); // important to change precondiconning for solving
//SOLVE THE PROBLEM
uh[] = A^-1* b; // solve the problem with more locking condition
//plot(uh, cmm = "uh after solver", wait=debug);
lambda[] = AA * uh[]; // compute the residual ( fast with matrix)
lambda[] = -rhs-lambda[]; // remark rhs -
Ik = (lambda + c*(uh-pay)) < 0; // set the new value
//plot(Ik, wait=debug,cmm=" lock set ",value=1 );
//plot(uh,wait=debug,cmm="uh");
// trick to compute $L^2$ norm of the variation
real[int] diff(n),Mdiff(n);
diff= uh[]-uhp[];
Mdiff = M*diff;
real err = sqrt(Mdiff'*diff);
cout << " ----- err norm L2 " << err << " ----- kkadapt = " << kkadapt <<endl;
if(err<eps && kkadapt )
{
break; //stop test
cout <<" STOP test at k = "<<k<<endl;
}
bool adapt = err<eps || (j>T/dt/4); //sure for each ...
if(adapt && iter>1)
{
kadapt++;
j=0;
Th=adaptmesh(Th,uh,err=tol);
kkadapt = tol == tolmin; // we reacht the bound
tol = max(tol/2,tolmin);
cout << " ++ tol = " << tol << " kadapt = " << kadapt << " reachted bound? " << kkadapt <<endl;
//plot(Th,wait=debug);
//plot(uh, wait=debug);
Vh xveloc = -x*r+x*sigmax^2+x*rho*sigmax*sigmay/2;
Vh yveloc = -y*r+y*sigmay^2+y*rho*sigmax*sigmay/2;
vh=max(uh,pay)-pay;
//plot(vh, wait=debug, cmm= "vh");
uhp=uhp;
n=uhp.n;
uh=uh;
Ik=Ik;
lambda=lambda;
pay=e;
A=a(Vh,Vh,tgv=tgv,solver=CG);
AA=a(Vh,Vh);
M=vM(Vh,Vh);
Aii.resize(n);
Aiin.resize(n);
Ak.resize(n);
b.resize(n);
rhs.resize(n);
Aii=A.diag; // get the diagonal of the matrix
}
k++;
cout<<"k = "<<k<<" iter = "<<iter<<endl;
}
j++;
iter++;
cout << "iter = "<< iter <<endl;
uhp[]=uh[] ; // set the previous solution
}
cout << endl;
cout << "===============================================" << endl;
cout << "==== CPU time =====" << endl;
cout << "===============================================" << endl;
cout << " ALL solving steps :::: " << clock()-time << endl;
cout<<endl;
cout << "==== uh(0,0)"<<uh(0,0)<<" =====" << endl;
cout<<endl;
cout << "==== uh(35,41)"<<uh(35,41)<<" =====" << endl;
vh=max(uh,pay)-pay;
plot(vh,wait=1,value=1,fill=1,cmm="vh");
plot(vh,bb=[[20,20],[50,50]],value=1, fill=0,cmm="vh-zoomed",wait=1,nbiso=50);
plot(uh,wait=1,value=1,fill=1,cmm="uh");
plot(uh,wait=1,value=1,fill=0,cmm="uh", nbiso=50);
plot(Th,cmm="Th",wait=1);
plot(Th,bb=[[20,20],[50,50]],cmm="Th-zoomed",wait=1);
savemesh(Th,"mm",[x,y,uh*10]);
// compute a cut
int i;
real[int] xx(L),yy(LL);
for (i=0;i<L;i++)
{
x=i/L; y=i/LL;
xx[i]=i;
yy[i]=uh(i,i); // value of uh at point (i/L , i/LL)
}
plot([xx,yy],bb=[[0,0],[80,80]], cmm="u like 1d", value = true, ps="VIA-american.eps",wait=true, aspectratio=true); // like gnuplot plot a cut of u
{ // File for gnuplot
ofstream gnu("plotVIA-american.gp");
for (int i = 0; i < L; i++)
gnu << xx[i] << " " << yy[i] << endl;
}
