Clear[τ,ρ,ξ,σ,r,q,t, T,S,K,η,M,ψ,Φ,i,j,n,nn,ii,A,η,η0,Φ];
(*some useful functions*)
Φ[ξ_]:=(1/2)+(1/2)Erf[(ξ/(Sqrt[2]))]; η[τ_,ρ_,T_]:=If[ρ==0,(τ/T),(1/( ρ T))(1-(E^(-ρ τ)))];η0[τ_,ρ_,T_]=If[ρ==0,((τ^3)/(3(T^2))),-((3+(E^(-2 ρ τ))-4 (E^(-ρ τ))-2 ρ τ)/(2 (T^2) (ρ^3)))]; 
nn=12;ii=6; (*this defines the number of expansion terms in α and τ*)
(*algorithm for caclulating α and τ expansion constants starts here*)
M=Table[0,{n,nn+1},{i,ii+1},{j,3(n-1)-1+4}];  
A[n_,i_,j_]:=If[j≤3 n,M[[1+n]][[1+i]][[5+j-2]],0]-(1+2 j)If[j≤3 n-2 && n>0, M[[1+n]][[1+i]][[5+j]],0]+(j+1) (j+2)If[j≤3 n-4, M[[1+n]][[1+i]][[5+j+2]],0];
Ψ1[i_]:=(((-1))^i)((2^(i+2))-2 )/(i+3)!;
Ψ2[i_]:=(((-1))^i)(397-29 (2^(9+i))+118 (3^(6+i))-9 (4^(7+i))+(5^(7+i))+(80-5 (2^(8+i))+4 (3^(6+i))) i+4 (i^2))/4/(i+7)!;
Ψ3[i_]:=(((-1))^i)(-15987+7971 (2^(10+i))+(2^(32+3 i))-523 (3^(12+i))+31 (2^(10+i)) (3^(12+i))+613 (4^(11+i))-197 (5^(11+i))-(7^(12+i))+2 (-1889+323 (2^(11+i))-80 (3^(11+i))+(2^(9+i)) (3^(11+i))+297 (4^(9+i))-3 (5^(11+i))) i+12 (-25+11 (2^(9+i))+(2^(17+2 i))-(3^(10+i))) (i^2)+8 (-1+(2^(7+i))) (i^3))/(i+12)!;

For[i=0,i≤ii,++i,
M=ReplacePart[M,{{1,1,3}->1}]];
For[n=1,n≤nn,++n,For[i=0,i≤ii,++i,For[j=3n-2,j≥0,--j,
M=ReplacePart[M,{1+n,1+i,5+j}->
-(3/((3(2+j-3 n) +2 (i+5 n) )))((∑_(k=1)^i)((2+j-3 n)(k+3)+2 (i-k+5 n) )Ψ1[k]M[[1+n]][[1+i-k]][[5+j]]-(∑_(k=0)^i)(Ψ1[k](k+3)(2+3 j+(j^2)) If[4+j≤3 n,M[[1+n]][[1+i-k]][[5+j+2]],0]-2 n Ψ2[k]A[n-1,i-k,j-1] +n (n-1)Ψ3[k]A[n-2,i-k,j-2] ))]]]];
(*Taylor expansion constants are stored in M and we are ready to write the solution*)
(*this is the solution to (Subscript[ψ, τ])=(1/2)(((ξ-(1/(ρ T))(1-(E^(-ρ τ)))))^2)(Subscript[ψ, ξ,ξ]), ψ[0]=(SuperPlus[ξ])*)
ψ[ξ_,τ_,ρ_,T_,η0_]:=ξ Φ[(ξ/(Sqrt[η0]))]+(1/(Sqrt[2 π]))(E^(-((ξ^2)/(2 η0)))) (∑_(n=0)^nn)(1/( n!((η0)^((3n-1)/2))(T^(3 n))(τ^(-5n))))(If[n==0,1,0]+(∑_(i=0)^ii)(∑_(j=0)^(3 n-2)) M[[1+n]][[1+i]][[5+j]]If[j==0,1,((η0)^(-(j/2)))(ξ^j)]If[i==0,1,(((ρ τ))^i)]);
(*the price*)
AsianFixedCall[i_,t_,T_,r_,q_,σ_,K_,S_]:=S (E^(-q t))ψ[((i-K T)/(S T))(E^(-(r-q)t))+η[t,r-q,T], (σ^2)t,((r-q)/(σ^2)), (σ^2)T,η0[(σ^2)t,((r-q)/(σ^2)),(σ^2)T]];
S=1.9;K=2.0;r=0.05;σ=0.5;T=1.0;t=T;q=0;
{1,NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.05;σ=0.5;T=1.0;t=T;q=0;
{2,NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.1;K=2.0;r=0.05;σ=0.5;T=1.0;t=T;q=0;
{3,NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.02;σ=0.1;T=1.0;t=T;q=0;
{4,NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.18;σ=0.3;T=1.0;t=T;q=0;
{5,NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.0125;σ=0.25;T=2.0;t=T;q=0;
{6,NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.05;σ=0.5;T=2.0;t=T;q=0;
{7,NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}

"low volatility test cases"
S=2.0;K=2.0;r=0.02;σ=0.1;T=1.0;t=T;q=0;
{4,NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.02;σ=0.05;T=1.0;t=T;q=0;
{4A,NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.02;σ=0.01;T=1.0;t=T;q=0;
{4B,NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.02;σ=0.005;T=1.0;t=T;q=0;
{4C,NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.02;σ=0.001;T=1.0;t=T;q=0;
{4D,NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
"test cases q>r"
S=1.9;K=2.0;r=0.05;q=0.1;σ=0.5;T=1.0;t=T;
{(SuperStar[1]),NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.05;q=0.1;σ=0.5;T=1.0;t=T;
{(SuperStar[2]),NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.1;K=2.0;r=0.05;q=0.1;σ=0.5;T=1.0;t=T;
{(SuperStar[3]),NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.02;q=0.04;σ=0.1;T=1.0;t=T;
{(SuperStar[4]),NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.18;q=0.36;σ=0.3;T=1.0;t=T;
{(SuperStar[5]),NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.0125;q=0.025;σ=0.25;T=2.0;t=T;
{(SuperStar[6]),NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.05;q=0.1;σ=0.5;T=2.0;t=T;
{(SuperStar[7]),NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
"low volatility and q>r"
S=2.0;K=2.0;r=0.02;q=0.04;σ=0.1;T=1.0;t=T;
{(SuperStar[4]),NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.02;q=0.04;σ=0.05;T=1.0;t=T;
{4(SuperStar[A]),NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.02;q=0.04;σ=0.01;T=1.0;t=T;
{4(SuperStar[B]),NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.02;q=0.04;σ=0.005;T=1.0;t=T;
{4(SuperStar[C]),NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.02;q=0.04;σ=0.001;T=1.0;t=T;
{4 (SuperStar[D]),NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
"test cases for q=r"
S=1.9;K=2.0;r=0.05;q=0.05;σ=0.5;T=1.0;t=T;
{"(1^( ** ))",NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.05;q=0.05;σ=0.5;T=1.0;t=T;
{"(2^( ** ))",NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.1;K=2.0;r=0.05;q=0.05;σ=0.5;T=1.0;t=T;
{"(3^( ** ))",NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.02;q=0.02;σ=0.1;T=1.0;t=T;
{"(4^( ** ))",NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.18;q=0.18;σ=0.3;T=1.0;t=T;
{"(5^( ** ))",NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.0125;q=0.0125;σ=0.25;T=2.0;t=T;
{"(6^( ** ))",NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.05;q=0.05;σ=0.5;T=2.0;t=T;
{"(7^( ** ))",NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}

"low volatility and q=r"
S=2.0;K=2.0;r=0.02;q=0.02;σ=0.1;T=1.0;t=T;
{"(4^( ** ))",NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.02;q=0.02;σ=0.05;T=1.0;t=T;
{"4(A^( ** ))",NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.02;q=0.02;σ=0.01;T=1.0;t=T;
{"4(B^( ** ))",NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.02;q=0.02;σ=0.005;T=1.0;t=T;
{"4(C^( ** ))",NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
S=2.0;K=2.0;r=0.02;q=0.02;σ=0.001;T=1.0;t=T;
{"4(D^( ** ))",NumberForm[AsianFixedCall[0,t,T,r,q,σ,K,S],10]}
Clear[τ,ρ,ξ,σ,r,q,t, T,S,K,η,M,ψ,Φ,i,j,n,nn,ii,A,η,η0,Φ];