Clear[c,b,B,Z,H,T,c1,nn,jn,j0,S,λ,r,σ,t,T,i,j,k,n,ξ,ζ,τ,η,Φ]
(*some useful functions*)
Φ[ξ_]:=(1/2)+(1/2)Erf[(ξ/(Sqrt[2]))]; η[τ_,ρ_,T_]:=If[ρ==0,(τ/T),(1/( ρ T))(1-(E^(-ρ τ)))];
(*define B,Z,H,T transformations, see [1] for details*)
B[c_]:=Module[{b=Table[0,{i,Length[c]},{j,Length[(c^\[Transpose])]},{l,j}]},
For[l=1,l<Length[(c^\[Transpose])],++l,For[j=l,j<Length[(c^\[Transpose])],++j,For[i=0,i<Length[c]-2j,++i,
b=ReplacePart[b,{1+i,1+j,1+j-l}->
(1/((i+j+1)))(If[i==0,0,b[[1+i-1]][[1+j]][[1+j-l]]]+((1/2)(i+1)(i+2) b[[1+i+2]][[1+j-1]][[1+j-l]]- (i+1) b[[1+i+1]][[1+j-1]][[1+j-l]]) +c[[1+i]][[1+j]][[1+j-l]])]]]];b=Drop[b,-2(j-1)]];
Z[c_]:=Module[{b=Table[0,{i,Length[c]},{j,Length[(c^\[Transpose])]},{k,j}]},
For[i=0,i<Length[b]-2,++i,For[j=0,j<Length[(b^\[Transpose])],++j,For[k=0,k<j,++k,
b=ReplacePart[b,{i+1,j+1,k+1}->If[i≥2,c[[1+i-2]][[1+j]][[1+k]],0]+If[j≥1&& k<j&&i≥1,c[[1+i-1]][[1+j-1]][[1+k]],0]-
If[j≥1&& k<j,(1+2 i)c[[1+i]][[1+j-1]][[1+k]],0]-
If[j≥2&&k<j-1, (i+1)c[[1+i+1]][[1+j-2]][[1+k]],0]+
If[j≥2&&k<j-1,(i+2) (i+1)c[[1+i+2]][[1+j-2]][[1+k]],0]];]]];b=Drop[b,-2]];
H[c_]:=Module[{b=Table[0,{i,Length[c]},{j,Length[(c^\[Transpose])]},{k,j}],c1=c},
For[i=0,i<Length[b],++i,For[j=0,j<Length[(b^\[Transpose])],++j,For[k=0,k<j,++k,
b=ReplacePart[b,{i+1,j+1,k+1}->(∑_(n=0)^k)((((-1))^(k-n))/((k-n+1)!))c1[[1+i]][[1+n+j-k]][[1+n]]];]]];c1=b;
For[i=0,i<Length[b],++i,For[j=0,j<Length[(b^\[Transpose])],++j,For[k=0,k<j,++k,
b=ReplacePart[b,{i+1,j+1,k+1}->(∑_(m=k+1)^j)((2^(-j+m))/((j-m)!))c1[[1+i]][[1+m]][[1+k]]];]]];c1=b;
For[i=0,i<Length[b],++i,For[j=0,j<Length[(b^\[Transpose])],++j,For[k=0,k<j,++k,
b=ReplacePart[b,{i+1,j+1,k+1}->(∑_(l=0)^i)((((-1))^(i-l))/((i-l)!))c1[[1+l]][[1+j]][[1+k]]];]]];b];
T[c_]:=Module[{b=Table[0,{i,Length[c]},{j,Length[(c^\[Transpose])]},{k,j}]},
For[i=0,i<Length[b],++i,For[j=1,j<Length[(b^\[Transpose])],++j,For[k=0,k<j,++k,
b=ReplacePart[b,{i+1,j+1,k+1}->c[[1+i]][[1+j-1]][[1+k]]];]]];b];
nn=7;jn=8;j0=jn;i0=2jn+2(jn+1)nn-2; (*this defines number of expansion terms in β, τ and ζ, see [1] for details*)
(*algorithm for caclulating β, τ and ζ expansion constants starts here*)
c1=Table[0,{i,i0+1},{j,j0+1},{k,j}];  For[j=1,j≤ j0,++j,c1=ReplacePart[c1,{1,1+j,1+j-1}-> ((((-1))^j)/(j!))]];
b={B[c1]}; If[nn>1,b=Join[b,{B[H[-2Z[b[[1]]]-T[Drop[c1,-(2jn+2)]]]]}]];
For[n=3,n≤nn,++n, b=Join[b,{B[H[-n Z[b[[n-1]]]+(1/2)n(n-1)T[H[Z[Drop[b[[n-2]],-(2jn+2)]]]]]]}];Print[StringJoin[ToString[n]," β orders calculated."]]];
For[n=nn,n≥1,--n, b[[n]]=Drop[b[[n]],-2(jn+1)(nn-n)]]
(*Taylor expansion constants are stored in b and we are ready to write the solution*)
(*this is the solution to (Subscript[ψ, τ])=(1/2)(((ξ-(1/(ρ T))(1-(E^(-ρ τ)))))^2)(Subscript[ψ, ξ,ξ]), ψ[0]=(SuperPlus[(ξ-λ)]) and ζ=Log[(ξ/λ)]+(1/2)τ, it is denoted (SubPlus[ζ]) in [1]*)

ψ[ζ_,τ_,ρ_,T_,λ_]:=(E^(ζ-(τ/2))) λ Φ[(ζ/(Sqrt[τ]))] -λ Φ [((ζ-τ)/( (Sqrt[τ])))]+(1/(Sqrt[2π τ]))((τ/T))(E^(-((ζ^2)/(2τ)))) (∑_(n=1)^nn)(∑_(i=0)^(Length[b[[1]]]-1))(∑_(j=0)^(Length[((b[[1]])^\[Transpose])]-1))(∑_(k=0)^(j-1))(1/(n!(((λ T))^(n-1))))b[[n]][[1+i]][[1+j]][[1+k]](ζ^i)(τ^j)If[k==0,1,(ρ^k)]

(*the price*)
AsianFloatingPut[i_,t_,T_,r_,q_,σ_,λ_,S_]:=S (E^(-q t)) ψ[Log[(1/λ)((i/(S T)) (E^(-(r-q) t))+η[t,r-q,T])]+(1/2) (σ^2) t ,(σ^2) t,(((r-q))/(σ^2)),(σ^2) T,λ];
"comparison of results"
S=1.9;K=2.0;r=0.05;σ=0.5;T=1.0;
{1,NumberForm[AsianFloatingPut[0,T,T,0,r,σ,(K/S),S],3]}
(*For[k=0,k≤25,++k,Print[NumberForm[AsianFloatingPut1[k,0,T,T,0,r,σ,(K/S),S],3]]]*)
S=2.0;K=2.0;r=0.05;σ=0.5;T=1.0;
{2,NumberForm[AsianFloatingPut[0,T,T,0,r,σ,(K/S),S],3]}
S=2.1;K=2.0;r=0.05;σ=0.5;T=1.0;
{3,NumberForm[AsianFloatingPut[0,T,T,0,r,σ,(K/S),S],3]}
S=2.0;K=2.0;r=0.02;σ=0.1;T=1.0;
{4,NumberForm[AsianFloatingPut[0,T,T,0,r,σ,(K/S),S],3]}
S=2.0;K=2.0;r=0.18;σ=0.3;T=1.0;
{5,NumberForm[AsianFloatingPut[0,T,T,0,r,σ,(K/S),S],3]}
S=2.0;K=2.0;r=0.0125;σ=0.25;T=2.0;
{6,NumberForm[AsianFloatingPut[0,T,T,0,r,σ,(K/S),S],3]}
S=2.0;K=2.0;r=0.05;σ=0.5;T=2.0;
{7,NumberForm[AsianFloatingPut[0,T,T,0,r,σ,(K/S),S],3]}

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

"low volatility and q=r"
S=2.0;K=2.0;r=0.02;q=0.02;σ=0.1;T=1.0;
{"(4^( ** ))",NumberForm[AsianFloatingPut[0,T,T,q,r,σ,(K/S),S],3]}
S=2.0;K=2.0;r=0.02;q=0.02;σ=0.05;T=1.0;
{"4(A^( ** ))",NumberForm[AsianFloatingPut[0,T,T,q,r,σ,(K/S),S],3]}
S=2.0;K=2.0;r=0.02;q=0.02;σ=0.01;T=1.0;
{"4(B^( ** ))",NumberForm[AsianFloatingPut[0,T,T,q,r,σ,(K/S),S],3]}
S=2.0;K=2.0;r=0.02;q=0.02;σ=0.005;T=1.0;
{"4(C^( ** ))",NumberForm[AsianFloatingPut[0,T,T,q,r,σ,(K/S),S],3]}
S=2.0;K=2.0;r=0.02;q=0.02;σ=0.001;T=1.0;
{"4(D^( ** ))",NumberForm[AsianFloatingPut[0,T,T,q,r,σ,(K/S),S],3]}
Clear[c,b,B,Z,H,T,c1,nn,jn,j0,S,λ,r,σ,t,T,i,j,k,n,ξ,ζ,τ,η,Φ];