Calculating the optimal portfolio numerically

Define the densities

qx[x_,sigma1_,t_] := PDF[NormalDistribution[sigma1^2*t,sigma1*Sqrt[t]],x]

qy[x_,sigma1_,t_] := PDF[NormalDistribution[0,sigma1*Sqrt[t]],x]

qm[x_,mu_,sigma2_,t_] := PDF[NormalDistribution[mu*t,sigma2*Sqrt[t]],x]

pmax1[x_,k1_,sigma1_,t_,x0_,y0_]:= Max[qx[x,sigma1,t]*x0, k1*qy[x,sigma1,t]*y0]/(x0*(1-CDF[NormalDistribution[sigma1^2*t,sigma1*Sqrt[t]], sigma1^2*t/2-Log[x0/(k1*y0)]])+k1*y0(CDF[NormalDistribution[0,sigma1*Sqrt[t]], sigma1^2*t/2-Log[x0/(k1*y0)]]))

No options portfolio

HardMaximum[mu_,sigma1_,sigma2_,t_]:= FindMaximum[{NIntegrate[qm[x,mu,sigma2,t]*Log[(w*qx[x,sigma1,t ]+ (1-w)*qy[x,sigma1,t])/qy[x,sigma1,t]],{x,-Infinity,Infinity}]},{w}]

HardMaximum[0.03,0.2,0.2,10]

HardMaximumPayoff[mu_,sigma1_,sigma2_,t_,w_]:= NIntegrate[qm[x,mu,sigma2,t]*Log[(w*qx[x,sigma1,t] + (1-w)qy[x,sigma1,t])/qy[x,sigma1,t]],{x,-Infinity,Infinity}]

One contract portfolio

HardMaximum1[mu_,sigma1_,sigma2_,t_,x0_,y0_]:= FindMaximum[{NIntegrate[qm[x,mu,sigma2,t]*Log[(w1*qx[x,sigma1,t ]+ w2*qy[x,sigma1,t] + (1-w1-w2)pmax1[x,k1,sigma1,t,x0,y0])/qy[x,sigma1,t]],{x,-Infinity,Infinity}],k1>0.1},{w1,w2,k1}]

HardMaximum1[0.05,0.3,0.3,10,1,1]

HardMaximum1[0.05,0.3,0.35,1,1,1]

2 contracts portfolio

HardMaximum2[mu_,sigma1_,sigma2_,t_,x0_,y0_]:= NMaximize[{Re[NIntegrate[qm[x,mu,sigma2,t]*Log[(w1*qx[x,sigma1,t ]+ w2*qy[x,sigma1,t] +  w3*pmax1[x,k1,sigma1,t,x0,y0]+ (1-w1-w2-w3)*pmax1[x,k2,sigma1,t,x0,y0])/qy[x,sigma1,t]],{x,-Infinity,Infinity}]],5>k1>0.1 && 5>k2>0.1 && w1+w2+w3 ==1 && Im[NIntegrate[qm[x,mu,sigma2,t]*Log[(w1*qx[x,sigma1,t ]+ w2*qy[x,sigma1,t] +  w3*pmax1[x,k1,sigma1,t,x0,y0]+ (1-w1-w2-w3)*pmax1[x,k2,sigma1,t,x0,y0])/qy[x,sigma1,t]],{x,-Infinity,Infinity}]] == 0},{w1,w2,w3,k1,k2}]

HardMaximum2[0.05,0.3,0.4,10,1,1]

HardPayoff2[mu_,sigma1_,sigma2_,t_,x0_,y0_, k1_,k2_,w1_,w2_,w3_]:= NIntegrate[qm[x,mu,sigma2,t]*Log[(w1*qx[x,sigma1,t ]+ w2*qy[x,sigma1,t] +  w3*pmax1[x,k1,sigma1,t,x0,y0]+ (1-w1-w2-w3)*pmax1[x,k2,sigma1,t,x0,y0])/qy[x,sigma1,t]],{x,-Infinity,Infinity}]

Equal contract weights portfolio

HardMaximum2SW[mu_,sigma1_,sigma2_,t_,x0_,y0_]:= NMaximize[{NIntegrate[qm[x,mu,sigma2,t]*Log[(w1*qx[x,sigma1,t ]+ (1-w1-2w2)*qy[x,sigma1,t] +  w2*pmax1[x,k1,sigma1,t,x0,y0]+ w2*pmax1[x,k2,sigma1,t,x0,y0])/qy[x,sigma1,t]],{x,-Infinity,Infinity}],k1>0.1 && k2>0.1&& w1+2*w2== 1},{w1,w2,k1,k2}]

Markowitz mean variance approach

Define densities

px[x_,sigma1_,t_] := PDF[NormalDistribution[sigma1^2*t,sigma1*Sqrt[t]],x]

py[x_,sigma1_,t_] := PDF[NormalDistribution[0,sigma1*Sqrt[t]],x]

pm[x_,mu_,sigma2_,t_] := PDF[NormalDistribution[mu*t,sigma2*Sqrt[t]],x]

State price density of an asset with payoff function  .

pmax[x_,k_,sigma1_,t_,x0_,y0_]:= Max[px[x,sigma1,t]*x0, k*py[x,sigma1,t]*y0]/(x0*(1-CDF[NormalDistribution[sigma1^2*t,sigma1*Sqrt[t]], sigma1^2*t/2-Log[x0/(k*y0)]])+k*y0(CDF[NormalDistribution[0,sigma1*Sqrt[t]], sigma1^2*t/2-Log[x0/(k*y0)]]))

Max contract density plot

Plot[pmax[x,2,0.2,20,1,1],{x,-5,5},PlotRange->Full]

Integrate[pmax[x,2,0.2,20,1,1],{x,-Infinity,Infinity}]

Entries of the  matrix:
(1,1)

Integrate[(px[x,sigma1,t]/py[x,sigma1,t])^2*py[x,sigma1,t],{x,-Infinity,Infinity}]

Simplify[((E^((sigma1^2) t))/(sigma1 (Sqrt[(1/((sigma1^2) t))]) (Sqrt[t]))),{sigma2>0,t>0,sigma1>0}]

p11[sigma1_,sigma2_,t_,mu_]:=(E^((sigma1^2) t))-1

(1,2)

Integrate[px[x,sigma1,t]*py[x,sigma1,t]/(py[x,sigma1,t])^2*py[x,sigma1,t],{x,-Infinity,Infinity}]

p12[sigma1_,sigma2_,t_,mu_] := 0

(2,2)

Integrate[py[x,sigma1,t]*py[x,sigma1,t]/(py[x,sigma1,t])^2*py[x,sigma1,t],{x,-Infinity,Infinity}]

p22[sigma1_,sigma2_,t_,mu_]:=0

(1,3)

px[x,sigma1,t]*pmax[x,k,sigma1,t,x0,y0]/(py[x,sigma1,t])^2*py[x,sigma1,t]

p1max[x_,k_,sigma1_,sigma2_,t_,x0_,y0_,mu_] := (((E^(((x^2)/(2 (sigma1^2) t))-((((-(sigma1^2) t+x))^2)/(2 (sigma1^2) t))))(((E^(-((x^2)/(2 (sigma1^2) t)))) k y0)/((Sqrt[2 π]) sigma1 (Sqrt[t]))))/( ((1/2) k y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))))

p2max[x_,k_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:= (((E^(((x^2)/(2 (sigma1^2) t))-((((-(sigma1^2) t+x))^2)/(2 (sigma1^2) t)))) (((E^(-((((-(sigma1^2) t+x))^2)/(2 (sigma1^2) t)))) x0)/((Sqrt[2 π]) sigma1 (Sqrt[t]))))/( ((1/2) k y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))))

Simplify[Integrate[p1max[x,k,sigma1,sigma2,t,x0,y0,mu],{x,-Infinity,(t*sigma1^2)/2-Log[x0/(y0*k)]}]+Integrate[p2max[x,k,sigma1,sigma2,t,x0,y0,mu],{x,(t*sigma1^2)/2-Log[x0/(y0*k)],Infinity}],{sigma1>0,t>0,sigma2>0}]

p13[k_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:= (((E^((sigma1^2) t)) x0+k y0-k y0 Erf[(((sigma1^2) t+2 Log[(x0/(k y0))])/(2 (Sqrt[2]) sigma1 (Sqrt[t])))]+(E^((sigma1^2) t)) x0 Erf[((3 (sigma1^2) t+2 Log[(x0/(k y0))])/(2 (Sqrt[2]) sigma1 (Sqrt[t])))])/(x0+x0 Erf[(((((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+k y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))-1

(2,3)

py[x,sigma1,t]*pmax[x,k,sigma1,t,x0,y0]/(py[x,sigma1,t])^2*py[x,sigma1,t]

p23[k_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:=0

(3,3)

(pmax[x,k,sigma1,t,x0,y0]/py[x,sigma1,t])^2*py[x,sigma1,t]

p1max[x_,k_,sigma1_,sigma2_,t_,x0_,y0_,mu_] := (((E^((x^2)/(2 (sigma1^2) t))) (Sqrt[2 π]) sigma1 (Sqrt[t])((((E^(-((x^2)/(2 (sigma1^2) t)))) k y0)/((Sqrt[2 π]) sigma1 (Sqrt[t]))))^2)/( ((1/2) k y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))^2))

p2max[x_,k_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:= (((E^((x^2)/(2 (sigma1^2) t))) (Sqrt[2 π]) sigma1 (Sqrt[t])((((E^(-((((-(sigma1^2) t+x))^2)/(2 (sigma1^2) t)))) x0)/((Sqrt[2 π]) sigma1 (Sqrt[t]))))^2)/(((1/2) k y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))^2))

Simplify[Integrate[p1max[x,k,sigma1,sigma2,t,x0,y0,mu],{x,-Infinity,(t*sigma1^2)/2-Log[x0/(y0*k)]}]+Integrate[p2max[x,k,sigma1,sigma2,t,x0,y0,mu],{x,(t*sigma1^2)/2-Log[x0/(y0*k)],Infinity}],{sigma1>0,t>0,sigma2>0}]

p33[k_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:=((2 ((E^((sigma1^2) t)) (x0^2)+(k^2) (y0^2)+(k^2) (y0^2) Erf[(((sigma1^2) t-2 Log[(x0/(k y0))])/(2 (Sqrt[2]) sigma1 (Sqrt[t])))]+(E^((sigma1^2) t)) (x0^2) Erf[((3 (sigma1^2) t+2 Log[(x0/(k y0))])/(2 (Sqrt[2]) sigma1 (Sqrt[t])))]))/(((k y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]-x0 (-2+Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])))^2))-1

(3,4),  WLOG :

Simplify[Integrate[((((E^(-((x^2)/(2 (sigma1^2) t)))) k1 y0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))/((1/2) k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])))*((((E^(-((x^2)/(2 (sigma1^2) t)))) k2 y0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))/((1/2) k2 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])))/py[x,sigma1,t],{x,-Infinity,(t*sigma1^2)/2-Log[x0/(y0*k1)]},Assumptions-> k1<k2]+Integrate[(((((E^(-((((-(sigma1^2) t+x))^2)/(2 (sigma1^2) t)))) x0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))/((1/2) k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))))*((((E^(-((x^2)/(2 (sigma1^2) t)))) k2 y0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))/((1/2) k2 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])))/py[x,sigma1,t],{x,(t*sigma1^2)/2-Log[x0/(y0*k1)],(t*sigma1^2)/2-Log[x0/(y0*k2)]},Assumptions-> k1<k2]+Integrate[(((((E^(-((((-(sigma1^2) t+x))^2)/(2 (sigma1^2) t)))) x0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))/((1/2) k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))))*(((((E^(-((((-(sigma1^2) t+x))^2)/(2 (sigma1^2) t)))) x0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))/((1/2) k2 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))))/py[x,sigma1,t],{x,(t*sigma1^2)/2-Log[x0/(y0*k2)],Infinity},Assumptions-> k1<k2],{sigma1>0,t>0,k1>0,x0>0,y0>0, k1<k2}]

p34[mu_,k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_]:=2 (((k2 x0 y0 (Erf[(((sigma1^2) t+2 Log[(x0/(k1 y0))])/(2 (Sqrt[2]) sigma1 (Sqrt[t])))]-Erf[(((sigma1^2) t+2 Log[(x0/(k2 y0))])/(2 (Sqrt[2]) sigma1 (Sqrt[t])))]))/((x0+x0 Erf[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]) (x0+x0 Erf[(((((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+k2 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])))+((k1 k2 (y0^2) (1+Erf[(((sigma1^2) t-2 Log[(x0/(k1 y0))])/(2 (Sqrt[2]) sigma1 (Sqrt[t])))]))/((k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]-x0 (-2+Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])) (k2 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]-x0 (-2+Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))))+(((E^((sigma1^2) t)) (x0^2) (1+Erf[((3 (sigma1^2) t+2 Log[(x0/(k2 y0))])/(2 (Sqrt[2]) sigma1 (Sqrt[t])))]))/((k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]-x0 (-2+Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])) (k2 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]-x0 (-2+Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])))))-1

DRIFT (under )

mu1[sigma1_,sigma2_,mu_,t_]:= (E^((1/2) (2 mu-(sigma1^2)+(sigma2^2)) t))-1

mu2[sigma1_,sigma2_,mu_,t_]:= 0

pmax[x,k,sigma1,t,x0,y0]*pm[x,mu,sigma2,t]/py[x,sigma1,t]

p1max[x_,k_,sigma1_,sigma2_,t_,x0_,y0_,mu_] := (((E^(((x^2)/(2 (sigma1^2) t))-((((-mu t+x))^2)/(2 (sigma2^2) t)))) sigma1 ((((E^(-((x^2)/(2 (sigma1^2) t)))) k y0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))))/( sigma2((1/2) k y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))))

p2max[x_,k_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:= (((E^(((x^2)/(2 (sigma1^2) t))-((((-mu t+x))^2)/(2 (sigma2^2) t)))) sigma1 ((((E^(-((((-(sigma1^2) t+x))^2)/(2 (sigma1^2) t)))) x0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))))/(sigma2((1/2) k y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))))

Simplify[Integrate[p1max[x,k,sigma1,sigma2,t,x0,y0,mu],{x,-Infinity,(t*sigma1^2)/2-Log[x0/(y0*k)]}]+Integrate[p2max[x,k,sigma1,sigma2,t,x0,y0,mu],{x,(t*sigma1^2)/2-Log[x0/(y0*k)],Infinity}],{sigma1>0,t>0,sigma2>0}]

mu3[k_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:=-((k y0 (-1+Erf[((2 mu t-(sigma1^2) t+2 Log[(x0/(k y0))])/(2 (Sqrt[2]) sigma2 (Sqrt[t])))]))/(k y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]-x0 (-2+Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])))-((E^((1/2) (2 mu-(sigma1^2)+(sigma2^2)) t)) x0 (1+Erf[(((mu+(sigma2^2)) (Sqrt[t]))/((Sqrt[2]) sigma2))]-(((mu+(sigma2^2)) Erf[((Sqrt[(((((mu+(sigma2^2)))^2) t)/(sigma2^2))])/(Sqrt[2]))])/(Sqrt[(((mu+(sigma2^2)))^2)]))+((Erf[((Sqrt[(((((2 mu-(sigma1^2)+2 (sigma2^2)) t+2 Log[(x0/(k y0))]))^2)/((sigma2^2) t))])/(2 (Sqrt[2])))] ((2 mu-(sigma1^2)+2 (sigma2^2)) t+2 Log[(x0/(k y0))]))/(Sqrt[((((2 mu-(sigma1^2)+2 (sigma2^2)) t+2 Log[(x0/(k y0))]))^2)]))))/(-k y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (-2+Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))-1

mu3[1,0.3,0.5,1,1,1,0.05]

No contract portfolio

B0[k_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:= {{p11[sigma1,sigma2,t,mu],p12[sigma1,sigma2,t,mu]},{p12[sigma1,sigma2,t,mu],p22[sigma1,sigma2,t,mu]}}

B0[k,sigma1,sigma2,t,x0,y0,mu]

Find the weights:

Solve[D[w1*mu1[sigma1,sigma2,mu,t]   -1/2*{w1,w2}.B0[k,sigma1,sigma2,t,x0,y0,mu].{{w1},{w2}} + l*(w1+w2-1),w1]==0 &&
D[w1*mu1[sigma1,sigma2,mu,t]   -1/2*{w1,w2}.B0[k,sigma1,sigma2,t,x0,y0,mu].{{w1},{w2}} + l*(w1+w2-1),w2]==0
&& D[w1*mu1[sigma1,sigma2,mu,t]    -1/2*{w1,w2}.B0[k,sigma1,sigma2,t,x0,y0,mu].{{w1},{w2}} + l*(w1+w2-1),l]==0, {w1,w2,l}]

vaha[mu_,sigma1_,sigma2_,t_]:=((-1+(E^((1/2) (2 mu-(sigma1^2)+(sigma2^2)) t)))/(-1+(E^((sigma1^2) t))))

Limit as :

Limit[vaha[mu,sigma1,sigma1,t], t->0, Direction->"FromAbove"]

Weights and payoff for

vaha[0.03,0.2,0.2,10]

MarkowitzPayoff0[mu_,sigma1_,sigma2_,t_]:=NIntegrate[pm[x,mu,sigma2,t]*Log[(vaha[mu,sigma1,sigma2,t]*px[x,sigma1,t] + (1-vaha[mu,sigma1,sigma2,t])*py[x,sigma1,t])/py[x,sigma1,t]],{x,-Infinity,Infinity}]

MarkowitzDensity0[x_,mu_,sigma1_,sigma2_,t_]:=vaha[mu,sigma1,sigma2,t]*px[x,sigma1,t] + (1-vaha[mu,sigma1,sigma2,t])*py[x,sigma1,t]

MarkowitzPayoff0[0.03,0.2,0.2,10]

One contract portfolio

Find the optimal weights:

A:= {{p11,p12,p13},{p12,p22,p23},{p13,p23,p33}}

B[k_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:= {{p11[sigma1,sigma2,t,mu],p12[sigma1,sigma2,t,mu],p13[k,sigma1,sigma2,t,x0,y0,mu]},{p12[sigma1,sigma2,t,mu],p22[sigma1,sigma2,t,mu],p23[k,sigma1,sigma2,t,x0,y0,mu]},{p13[k,sigma1,sigma2,t,x0,y0,mu],p23[k,sigma1,sigma2,t,x0,y0,mu],p33[k,sigma1,sigma2,t,x0,y0,mu]}}

Solve[D[w1*mu1  + w3*mu3  -1/2*{w1,w2,w3}.A.{{w1},{w2},{w3}} + l*(w1+w2+w3-1),w1]==0&& D[w1*mu1  + w3*mu3  -1/2*{w1,w2,w3}.A.{{w1},{w2},{w3}} + l*(w1+w2+w3-1),w2]==0
&&D[w1*mu1  + w3*mu3  -1/2*{w1,w2,w3}.A.{{w1},{w2},{w3}} + l*(w1+w2+w3-1),w3]==0
&& D[w1*mu1  + w3*mu3  -1/2*{w1,w2,w3}.A.{{w1},{w2},{w3}} + l*(w1+w2+w3-1),l]==0, {w1,w2,w3,l}]

w1[k_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:= -((-((p12[sigma1,sigma2,t,mu]-p13[k,sigma1,sigma2,t,x0,y0,mu]-p22[sigma1,sigma2,t,mu]+p23[k,sigma1,sigma2,t,x0,y0,mu]) (-mu1[sigma1,sigma2,mu,t]+mu3[k,sigma1,sigma2,t,x0,y0,mu]+p13[k,sigma1,sigma2,t,x0,y0,mu]-p33[k,sigma1,sigma2,t,x0,y0,mu]))+(-mu1[sigma1,sigma2,mu,t]+p13[k,sigma1,sigma2,t,x0,y0,mu]-p23[k,sigma1,sigma2,t,x0,y0,mu]) (p12[sigma1,sigma2,t,mu]-p13[k,sigma1,sigma2,t,x0,y0,mu]-p23[k,sigma1,sigma2,t,x0,y0,mu]+p33[k,sigma1,sigma2,t,x0,y0,mu]))/(-((p12[sigma1,sigma2,t,mu]-p13[k,sigma1,sigma2,t,x0,y0,mu]-p22[sigma1,sigma2,t,mu]+p23[k,sigma1,sigma2,t,x0,y0,mu]) (p11[sigma1,sigma2,t,mu]-2 p13[k,sigma1,sigma2,t,x0,y0,mu]+p33[k,sigma1,sigma2,t,x0,y0,mu]))+(p11[sigma1,sigma2,t,mu]-p12[sigma1,sigma2,t,mu]-p13[k,sigma1,sigma2,t,x0,y0,mu]+p23[k,sigma1,sigma2,t,x0,y0,mu]) (p12[sigma1,sigma2,t,mu]-p13[k,sigma1,sigma2,t,x0,y0,mu]-p23[k,sigma1,sigma2,t,x0,y0,mu]+p33[k,sigma1,sigma2,t,x0,y0,mu])))

w2[k_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:=-((-mu3[k,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu]-mu1[sigma1,sigma2,mu,t] p12[sigma1,sigma2,t,mu]+mu3[k,sigma1,sigma2,t,x0,y0,mu] p12[sigma1,sigma2,t,mu]+mu1[sigma1,sigma2,mu,t] p13[k,sigma1,sigma2,t,x0,y0,mu]+mu3[k,sigma1,sigma2,t,x0,y0,mu] p13[k,sigma1,sigma2,t,x0,y0,mu]+p12[sigma1,sigma2,t,mu] p13[k,sigma1,sigma2,t,x0,y0,mu]-((p13[k,sigma1,sigma2,t,x0,y0,mu])^2)+mu1[sigma1,sigma2,mu,t] p23[k,sigma1,sigma2,t,x0,y0,mu]-mu3[k,sigma1,sigma2,t,x0,y0,mu] p23[k,sigma1,sigma2,t,x0,y0,mu]-p11[sigma1,sigma2,t,mu] p23[k,sigma1,sigma2,t,x0,y0,mu]+p13[k,sigma1,sigma2,t,x0,y0,mu] p23[k,sigma1,sigma2,t,x0,y0,mu]-mu1[sigma1,sigma2,mu,t] p33[k,sigma1,sigma2,t,x0,y0,mu]+p11[sigma1,sigma2,t,mu] p33[k,sigma1,sigma2,t,x0,y0,mu]-p12[sigma1,sigma2,t,mu] p33[k,sigma1,sigma2,t,x0,y0,mu])/(((p12[sigma1,sigma2,t,mu])^2)-2 p12[sigma1,sigma2,t,mu] p13[k,sigma1,sigma2,t,x0,y0,mu]+((p13[k,sigma1,sigma2,t,x0,y0,mu])^2)-p11[sigma1,sigma2,t,mu] p22[sigma1,sigma2,t,mu]+2 p13[k,sigma1,sigma2,t,x0,y0,mu] p22[sigma1,sigma2,t,mu]+2 p11[sigma1,sigma2,t,mu] p23[k,sigma1,sigma2,t,x0,y0,mu]-2 p12[sigma1,sigma2,t,mu] p23[k,sigma1,sigma2,t,x0,y0,mu]-2 p13[k,sigma1,sigma2,t,x0,y0,mu] p23[k,sigma1,sigma2,t,x0,y0,mu]+((p23[k,sigma1,sigma2,t,x0,y0,mu])^2)-p11[sigma1,sigma2,t,mu] p33[k,sigma1,sigma2,t,x0,y0,mu]+2 p12[sigma1,sigma2,t,mu] p33[k,sigma1,sigma2,t,x0,y0,mu]-p22[sigma1,sigma2,t,mu] p33[k,sigma1,sigma2,t,x0,y0,mu]))

Calculate optimal strike and payoff:

drift[k_,mu_,sigma1_,sigma2_,t_,x0_,y0_]:= w1[k,sigma1,sigma2,t,x0,y0,mu]*mu1[sigma1,sigma2,mu,t]  + (1-w1[k,sigma1,sigma2,t,x0,y0,mu]-w2[k,sigma1,sigma2,t,x0,y0,mu])*mu3[k,sigma1,sigma2,t,x0,y0,mu]

forma[k_,mu_,sigma1_,sigma2_,t_,x0_,y0_]:=-1/2*Norm[{w1[k,sigma1,sigma2,t,x0,y0,mu],w2[k,sigma1,sigma2,t,x0,y0,mu],1-w1[k,sigma1,sigma2,t,x0,y0,mu]-w2[k,sigma1,sigma2,t,x0,y0,mu]}.B[k,sigma1,sigma2,t,x0,y0,mu].{{w1[k,sigma1,sigma2,t,x0,y0,mu]},{w2[k,sigma1,sigma2,t,x0,y0,mu]},{1-w1[k,sigma1,sigma2,t,x0,y0,mu]-w2[k,sigma1,sigma2,t,x0,y0,mu]}}]

vynos[k_,mu_,sigma1_,sigma2_,t_,x0_,y0_]:= drift[k,mu,sigma1,sigma2,t,x0,y0] +forma[k,mu,sigma1,sigma2,t,x0,y0]

MarkowitzApMax1[mu_,sigma1_,sigma2_,t_,x0_,y0_]:= FindMaximum[{vynos[k,mu,sigma1,sigma2,t,x0,y0],k>0.1},{k}]

MarkowitzApMax1[0.05,0.3,0.3,10,1,1]

MarkowitzApMax1[0.05,0.3,0.35,1,1,1]

MarkowitzPayoff1[k_,mu_,sigma1_,sigma2_,t_,x0_,y0_]:=NIntegrate[pm[x,mu,sigma2,t]*Log[(w1[k,sigma1,sigma2,t,x0,y0,mu]*px[x,sigma1,t] + w2[k,sigma1,sigma2,t,x0,y0,mu]py[x,sigma1,t] + (1-w1[k,sigma1,sigma2,t,x0,y0,mu]-w2[k,sigma1,sigma2,t,x0,y0,mu])*pmax[x,k,sigma1,t,x0,y0])/py[x,sigma1,t]],{x,-Infinity,Infinity}]

MarkowitzDensity1[x_,k_,mu_,sigma1_,sigma2_,t_,x0_,y0_]:= (w1[k,sigma1,sigma2,t,x0,y0,mu]*px[x,sigma1,t] + w2[k,sigma1,sigma2,t,x0,y0,mu]py[x,sigma1,t] + (1-w1[k,sigma1,sigma2,t,x0,y0,mu]-w2[k,sigma1,sigma2,t,x0,y0,mu])*pmax[x,k,sigma1,t,x0,y0])

MarkowitzPayoff1[1.681,0.05,0.3,0.3,10,1,1]

Two contracts portfolio

Find the optimal weights:

A:= {{p11,p12,p13,p13k},{p12,p22,p23,p23k},{p13,p23,p33,p34},{p13k,p23k,p34,p33k}}

Solve[D[w1*mu1  + w3*mu3 + w4*mu3k  -1/2*{w1,w2,w3,w4}.A.{{w1},{w2},{w3},{w4}} + l*(w1+w2+w3+w4-1),w1]==0&& D[w1*mu1  + w3*mu3 + w4*mu3k  -1/2*{w1,w2,w3,w4}.A.{{w1},{w2},{w3},{w4}} + l*(w1+w2+w3+w4-1),w2]==0
&&D[w1*mu1  + w3*mu3 + w4*mu3k  -1/2*{w1,w2,w3,w4}.A.{{w1},{w2},{w3},{w4}} + l*(w1+w2+w3+w4-1),w3]==0 &&
D[w1*mu1  + w3*mu3 + w4*mu3k  -1/2*{w1,w2,w3,w4}.A.{{w1},{w2},{w3},{w4}} + l*(w1+w2+w3+w4-1),w4]==0
&& D[w1*mu1  + w3*mu3 + w4*mu3k  -1/2*{w1,w2,w3,w4}.A.{{w1},{w2},{w3},{w4}} + l*(w1+w2+w3+w4-1),l]==0, {w1,w2,w3,w4,l}]

Mw1[k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:=-(((-((p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]-p23[k1,sigma1,sigma2,t,x0,y0,mu]+p23[k2,sigma1,sigma2,t,x0,y0,mu]) (p11[sigma1,sigma2,t,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]-p23[k2,sigma1,sigma2,t,x0,y0,mu]+p33[k2,sigma1,sigma2,t,x0,y0,mu]))+(p11[sigma1,sigma2,t,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]-p11[sigma1,sigma2,t,mu]+p23[k2,sigma1,sigma2,t,x0,y0,mu]) (p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]+p33[k2,sigma1,sigma2,t,x0,y0,mu]-p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])) ((-mu1[sigma1,sigma2,t,mu]+mu3[k1,sigma1,sigma2,t,x0,y0,mu]+p13[k2,sigma1,sigma2,t,x0,y0,mu]-p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]) (p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]+p33[k2,sigma1,sigma2,t,x0,y0,mu]-p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])-(-mu1[sigma1,sigma2,t,mu]+mu3[k2,sigma1,sigma2,t,x0,y0,mu]+p13[k2,sigma1,sigma2,t,x0,y0,mu]-p33[k2,sigma1,sigma2,t,x0,y0,mu]) (p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]-p33[k1,sigma1,sigma2,t,x0,y0,mu]+p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]))-(-((p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]-p23[k1,sigma1,sigma2,t,x0,y0,mu]+p23[k2,sigma1,sigma2,t,x0,y0,mu]) (-mu1[sigma1,sigma2,t,mu]+mu3[k2,sigma1,sigma2,t,x0,y0,mu]+p13[k2,sigma1,sigma2,t,x0,y0,mu]-p33[k2,sigma1,sigma2,t,x0,y0,mu]))+(-mu1[sigma1,sigma2,t,mu]+p13[k2,sigma1,sigma2,t,x0,y0,mu]-p23[k2,sigma1,sigma2,t,x0,y0,mu]) (p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]+p33[k2,sigma1,sigma2,t,x0,y0,mu]-p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])) ((p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]+p33[k2,sigma1,sigma2,t,x0,y0,mu]-p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]) (p11[sigma1,sigma2,t,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]-p23[k1,sigma1,sigma2,t,x0,y0,mu]+p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])-(p11[sigma1,sigma2,t,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]-p23[k2,sigma1,sigma2,t,x0,y0,mu]+p33[k2,sigma1,sigma2,t,x0,y0,mu]) (p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]-p33[k1,sigma1,sigma2,t,x0,y0,mu]+p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])))/((-((p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]-p23[k1,sigma1,sigma2,t,x0,y0,mu]+p23[k2,sigma1,sigma2,t,x0,y0,mu]) (p11[sigma1,sigma2,t,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]-p23[k2,sigma1,sigma2,t,x0,y0,mu]+p33[k2,sigma1,sigma2,t,x0,y0,mu]))+(p11[sigma1,sigma2,t,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]-p11[sigma1,sigma2,t,mu]+p23[k2,sigma1,sigma2,t,x0,y0,mu]) (p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]+p33[k2,sigma1,sigma2,t,x0,y0,mu]-p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])) ((p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]+p33[k2,sigma1,sigma2,t,x0,y0,mu]-p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]) (p11[sigma1,sigma2,t,mu]-p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]+p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])-(p11[sigma1,sigma2,t,mu]-2 p13[k2,sigma1,sigma2,t,x0,y0,mu]+p33[k2,sigma1,sigma2,t,x0,y0,mu]) (p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]-p33[k1,sigma1,sigma2,t,x0,y0,mu]+p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]))-(-((p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]-p23[k1,sigma1,sigma2,t,x0,y0,mu]+p23[k2,sigma1,sigma2,t,x0,y0,mu]) (p11[sigma1,sigma2,t,mu]-2 p13[k2,sigma1,sigma2,t,x0,y0,mu]+p33[k2,sigma1,sigma2,t,x0,y0,mu]))+(p11[sigma1,sigma2,t,mu]-p11[sigma1,sigma2,t,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]+p23[k2,sigma1,sigma2,t,x0,y0,mu]) (p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]+p33[k2,sigma1,sigma2,t,x0,y0,mu]-p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])) ((p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]+p33[k2,sigma1,sigma2,t,x0,y0,mu]-p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]) (p11[sigma1,sigma2,t,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]-p23[k1,sigma1,sigma2,t,x0,y0,mu]+p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])-(p11[sigma1,sigma2,t,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]-p23[k2,sigma1,sigma2,t,x0,y0,mu]+p33[k2,sigma1,sigma2,t,x0,y0,mu]) (p13[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu]-p33[k1,sigma1,sigma2,t,x0,y0,mu]+p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]))))

Mw1[2,2,0.3,0.4,10,1,1,0.05]

Mw2[k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:= -((-mu3[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu]+mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu]-mu3[k2,sigma1,sigma2,t,x0,y0,mu] ((p13[k1,sigma1,sigma2,t,x0,y0,mu])^2)+mu3[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu]-mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu]+mu3[k1,sigma1,sigma2,t,x0,y0,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu]+mu3[k2,sigma1,sigma2,t,x0,y0,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu]-mu3[k1,sigma1,sigma2,t,x0,y0,mu] ((p13[k2,sigma1,sigma2,t,x0,y0,mu])^2)+mu3[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]-mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]-mu1[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]+mu3[k2,sigma1,sigma2,t,x0,y0,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]+mu1[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]-2 mu3[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]+mu3[k2,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]+p13[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]-((p13[k2,sigma1,sigma2,t,x0,y0,mu])^2) p23[k1,sigma1,sigma2,t,x0,y0,mu]-mu3[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+mu1[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+mu3[k1,sigma1,sigma2,t,x0,y0,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]-2 mu3[k2,sigma1,sigma2,t,x0,y0,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]-((p13[k1,sigma1,sigma2,t,x0,y0,mu])^2) p23[k2,sigma1,sigma2,t,x0,y0,mu]-mu1[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+mu3[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+p13[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]+mu1[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]-mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]-mu1[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]-mu3[k2,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]-p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]+((p13[k2,sigma1,sigma2,t,x0,y0,mu])^2) p33[k1,sigma1,sigma2,t,x0,y0,mu]-mu1[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]+mu3[k2,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]+p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]+mu3[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+mu1[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-mu3[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-mu1[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-mu3[k1,sigma1,sigma2,t,x0,y0,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-p11[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+((p13[k1,sigma1,sigma2,t,x0,y0,mu])^2) p33[k2,sigma1,sigma2,t,x0,y0,mu]-mu1[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+mu3[k1,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+mu1[sigma1,sigma2,t,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-p11[sigma1,sigma2,t,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+p11[sigma1,sigma2,t,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-mu3[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 mu1[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+mu3[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+mu1[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+mu3[k2,sigma1,sigma2,t,x0,y0,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+p11[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+mu1[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+mu3[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+mu1[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-mu3[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+mu1[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-mu3[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-mu1[sigma1,sigma2,t,mu] ((p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2)+p11[sigma1,sigma2,t,mu] ((p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2)-p11[sigma1,sigma2,t,mu] ((p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2))/(-((p13[k1,sigma1,sigma2,t,x0,y0,mu])^2) p11[sigma1,sigma2,t,mu]+2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu]-((p13[k2,sigma1,sigma2,t,x0,y0,mu])^2) p11[sigma1,sigma2,t,mu]+2 p11[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]-2 p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]-2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]+2 ((p13[k2,sigma1,sigma2,t,x0,y0,mu])^2) p23[k1,sigma1,sigma2,t,x0,y0,mu]-p11[sigma1,sigma2,t,mu] ((p23[k1,sigma1,sigma2,t,x0,y0,mu])^2)+2 p13[k2,sigma1,sigma2,t,x0,y0,mu] ((p23[k1,sigma1,sigma2,t,x0,y0,mu])^2)-2 p11[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+2 ((p13[k1,sigma1,sigma2,t,x0,y0,mu])^2) p23[k2,sigma1,sigma2,t,x0,y0,mu]+2 p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]-2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+2 p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]-2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]-2 p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]-p11[sigma1,sigma2,t,mu] ((p23[k2,sigma1,sigma2,t,x0,y0,mu])^2)+2 p13[k1,sigma1,sigma2,t,x0,y0,mu] ((p23[k2,sigma1,sigma2,t,x0,y0,mu])^2)-((p11[sigma1,sigma2,t,mu])^2) p33[k1,sigma1,sigma2,t,x0,y0,mu]+2 p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]-((p13[k2,sigma1,sigma2,t,x0,y0,mu])^2) p33[k1,sigma1,sigma2,t,x0,y0,mu]+p11[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]-2 p13[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]-2 p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]+2 p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]+2 p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]-((p23[k2,sigma1,sigma2,t,x0,y0,mu])^2) p33[k1,sigma1,sigma2,t,x0,y0,mu]-((p11[sigma1,sigma2,t,mu])^2) p33[k2,sigma1,sigma2,t,x0,y0,mu]+2 p11[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-((p13[k1,sigma1,sigma2,t,x0,y0,mu])^2) p33[k2,sigma1,sigma2,t,x0,y0,mu]+p11[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-2 p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+2 p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-((p23[k1,sigma1,sigma2,t,x0,y0,mu])^2) p33[k2,sigma1,sigma2,t,x0,y0,mu]+p11[sigma1,sigma2,t,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-2 p11[sigma1,sigma2,t,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+p11[sigma1,sigma2,t,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+2 ((p11[sigma1,sigma2,t,mu])^2) p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 p11[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 p11[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 p13[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 p23[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-p11[sigma1,sigma2,t,mu] ((p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2)+2 p11[sigma1,sigma2,t,mu] ((p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2)-p11[sigma1,sigma2,t,mu] ((p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2)))

Mw2[2,2,0.3,0.4,10,1,1,0.05]

Mw3[k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:=-((-mu3[k1,sigma1,sigma2,t,x0,y0,mu] ((p11[sigma1,sigma2,t,mu])^2)+mu3[k2,sigma1,sigma2,t,x0,y0,mu] ((p11[sigma1,sigma2,t,mu])^2)-mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu]+2 mu3[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu]-mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu]+mu3[k2,sigma1,sigma2,t,x0,y0,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu]-mu3[k1,sigma1,sigma2,t,x0,y0,mu] ((p13[k2,sigma1,sigma2,t,x0,y0,mu])^2)+mu3[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu]-mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu]-mu1[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu]+mu3[k2,sigma1,sigma2,t,x0,y0,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu]+mu1[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu]-2 mu3[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu]+mu3[k2,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu]+p13[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu]-((p13[k2,sigma1,sigma2,t,x0,y0,mu])^2) p11[sigma1,sigma2,t,mu]+mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]+mu1[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]-mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]-mu1[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]-mu3[k2,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]-p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]+((p13[k2,sigma1,sigma2,t,x0,y0,mu])^2) p23[k1,sigma1,sigma2,t,x0,y0,mu]-2 mu3[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]-mu1[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+2 mu3[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]-mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+2 mu1[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]-mu3[k2,sigma1,sigma2,t,x0,y0,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]-p11[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]-mu1[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+2 mu3[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+2 p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]-p13[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]-mu1[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+mu3[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]-p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+mu1[sigma1,sigma2,t,mu] ((p23[k2,sigma1,sigma2,t,x0,y0,mu])^2)-mu3[k1,sigma1,sigma2,t,x0,y0,mu] ((p23[k2,sigma1,sigma2,t,x0,y0,mu])^2)-p11[sigma1,sigma2,t,mu] ((p23[k2,sigma1,sigma2,t,x0,y0,mu])^2)+p13[k1,sigma1,sigma2,t,x0,y0,mu] ((p23[k2,sigma1,sigma2,t,x0,y0,mu])^2)+mu3[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+mu1[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-2 mu3[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-((p11[sigma1,sigma2,t,mu])^2) p33[k2,sigma1,sigma2,t,x0,y0,mu]-mu1[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+p11[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-mu1[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+mu3[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+p11[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-p13[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+mu1[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-mu1[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+((p11[sigma1,sigma2,t,mu])^2) p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+mu1[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+mu1[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-mu3[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-p11[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+p13[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-mu1[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])/(((p13[k1,sigma1,sigma2,t,x0,y0,mu])^2) p11[sigma1,sigma2,t,mu]-2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu]+((p13[k2,sigma1,sigma2,t,x0,y0,mu])^2) p11[sigma1,sigma2,t,mu]-2 p11[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]+2 p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]+2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu]-2 ((p13[k2,sigma1,sigma2,t,x0,y0,mu])^2) p23[k1,sigma1,sigma2,t,x0,y0,mu]+p11[sigma1,sigma2,t,mu] ((p23[k1,sigma1,sigma2,t,x0,y0,mu])^2)-2 p13[k2,sigma1,sigma2,t,x0,y0,mu] ((p23[k1,sigma1,sigma2,t,x0,y0,mu])^2)+2 p11[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]-2 ((p13[k1,sigma1,sigma2,t,x0,y0,mu])^2) p23[k2,sigma1,sigma2,t,x0,y0,mu]-2 p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]-2 p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+2 p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu]+p11[sigma1,sigma2,t,mu] ((p23[k2,sigma1,sigma2,t,x0,y0,mu])^2)-2 p13[k1,sigma1,sigma2,t,x0,y0,mu] ((p23[k2,sigma1,sigma2,t,x0,y0,mu])^2)+((p11[sigma1,sigma2,t,mu])^2) p33[k1,sigma1,sigma2,t,x0,y0,mu]-2 p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]+((p13[k2,sigma1,sigma2,t,x0,y0,mu])^2) p33[k1,sigma1,sigma2,t,x0,y0,mu]-p11[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]+2 p13[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]+2 p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]-2 p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]-2 p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu]+((p23[k2,sigma1,sigma2,t,x0,y0,mu])^2) p33[k1,sigma1,sigma2,t,x0,y0,mu]+((p11[sigma1,sigma2,t,mu])^2) p33[k2,sigma1,sigma2,t,x0,y0,mu]-2 p11[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+((p13[k1,sigma1,sigma2,t,x0,y0,mu])^2) p33[k2,sigma1,sigma2,t,x0,y0,mu]-p11[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+2 p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-2 p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+((p23[k1,sigma1,sigma2,t,x0,y0,mu])^2) p33[k2,sigma1,sigma2,t,x0,y0,mu]-p11[sigma1,sigma2,t,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]+2 p11[sigma1,sigma2,t,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-p11[sigma1,sigma2,t,mu] p33[k1,sigma1,sigma2,t,x0,y0,mu] p33[k2,sigma1,sigma2,t,x0,y0,mu]-2 ((p11[sigma1,sigma2,t,mu])^2) p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 p11[sigma1,sigma2,t,mu] p13[k1,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 p11[sigma1,sigma2,t,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p13[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 p11[sigma1,sigma2,t,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 p13[k2,sigma1,sigma2,t,x0,y0,mu] p11[sigma1,sigma2,t,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 p11[sigma1,sigma2,t,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 p13[k2,sigma1,sigma2,t,x0,y0,mu] p23[k1,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 p11[sigma1,sigma2,t,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 p13[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 p23[k1,sigma1,sigma2,t,x0,y0,mu] p23[k2,sigma1,sigma2,t,x0,y0,mu] p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+p11[sigma1,sigma2,t,mu] ((p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2)-2 p11[sigma1,sigma2,t,mu] ((p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2)+p11[sigma1,sigma2,t,mu] ((p34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2)))

Mw4[k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:= 1-Mw1[k1,k2,sigma1,sigma2,t,x0,y0,mu]-Mw2[k1,k2,sigma1,sigma2,t,x0,y0,mu]-Mw3[k1,k2,sigma1,sigma2,t,x0,y0,mu]

Mw4[2,3,0.3,0.3,10,1,1,0.05]

Calculate the optimal strikes and payoff:

MarkowitzPayoff2[k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:= NIntegrate[pm[x,mu,sigma2,t]*Log[(Mw1[k1,k2,sigma1,sigma2,t,x0,y0,mu]*px[x,sigma1,t] + Mw2[k1,k2,sigma1,sigma2,t,x0,y0,mu]*py[x,sigma1,t] + Mw3[k1,k2,sigma1,sigma2,t,x0,y0,mu]*pmax[x,k1,sigma1,t,x0,y0] +Mw4[k1,k2,sigma1,sigma2,t,x0,y0,mu]*pmax[x,k2,sigma1,t,x0,y0] )/py[x,sigma1,t]],{x,-Infinity,Infinity}]

MarkowitzDensity2[x_,k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:= (Mw1[k1,k2,sigma1,sigma2,t,x0,y0,mu]*px[x,sigma1,t] + Mw2[k1,k2,sigma1,sigma2,t,x0,y0,mu]*py[x,sigma1,t] + Mw3[k1,k2,sigma1,sigma2,t,x0,y0,mu]*pmax[x,k1,sigma1,t,x0,y0] +Mw4[k1,k2,sigma1,sigma2,t,x0,y0,mu]*pmax[x,k2,sigma1,t,x0,y0] )

B2[k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:= {{p11[sigma1,sigma2,t,mu],p12[sigma1,sigma2,t,mu],p13[k1,sigma1,sigma2,t,x0,y0,mu],p13[k2,sigma1,sigma2,t,x0,y0,mu]},{p12[sigma1,sigma2,t,mu],p22[sigma1,sigma2,t,mu],p23[k1,sigma1,sigma2,t,x0,y0,mu],p23[k2,sigma1,sigma2,t,x0,y0,mu]},{p13[k1,sigma1,sigma2,t,x0,y0,mu],p23[k1,sigma1,sigma2,t,x0,y0,mu],p33[k1,sigma1,sigma2,t,x0,y0,mu],p34[mu,k1,k2,sigma1,sigma2,t,x0,y0]},{p13[k2,sigma1,sigma2,t,x0,y0,mu],p23[k2,sigma1,sigma2,t,x0,y0,mu],p34[mu,k1,k2,sigma1,sigma2,t,x0,y0],p33[k2,sigma1,sigma2,t,x0,y0,mu]}}

drift2[k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:= mu1[sigma1,sigma2,t,mu]*Mw1[k1,k2,sigma1,sigma2,t,x0,y0,mu]  + Mw3[k1,k2,sigma1,sigma2,t,x0,y0,mu]*mu3[k1,sigma1,sigma2,t,x0,y0,mu]+ Mw4[k1,k2,sigma1,sigma2,t,x0,y0,mu]*mu3[k2,sigma1,sigma2,t,x0,y0,mu]

forma2[k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:= -1/2*Norm[{Mw1[k1,k2,sigma1,sigma2,t,x0,y0,mu],Mw2[k1,k2,sigma1,sigma2,t,x0,y0,mu],Mw3[k1,k2,sigma1,sigma2,t,x0,y0,mu],Mw4[k1,k2,sigma1,sigma2,t,x0,y0,mu]}.B2[k1,k2,sigma1,sigma2,t,x0,y0,mu].{{Mw1[k1,k2,sigma1,sigma2,t,x0,y0,mu]},{Mw2[k1,k2,sigma1,sigma2,t,x0,y0,mu]},{Mw3[k1,k2,sigma1,sigma2,t,x0,y0,mu]},{Mw4[k1,k2,sigma1,sigma2,t,x0,y0,mu]}}]

vynos2[k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:= drift2[k1,k2,sigma1,sigma2,t,x0,y0,mu]+forma2[k1,k2,sigma1,sigma2,t,x0,y0,mu]

MarkowitzApMax2[sigma1_,sigma2_,t_,x0_,y0_,mu_]:= FindMaximum[{vynos2[k1,k2,sigma1,sigma2,t,x0,y0,mu],20> k1>0.1 &&20>k2>0.1},{k1,k2}]

MarkowitzApMax2[0.3,0.4,10,1,1,0.05]

MarkowitzPayoff2[0.1,1.38792,0.3,0.4,10,1,1,0.05]

Mw1[0.1,1.38792,0.3,0.4,10,1,1,0.05]

Mw2[0.1,1.38792,0.3,0.4,10,1,1,0.05]

Mw3[0.1,1.38792,0.3,0.4,10,1,1,0.05]

MarkowitzPayoff2[0.1,1.38792,0.3,0.4,10,1,1,0.05]

Tuned mean-variance approach

Define the densities:

qx[x_,sigma1_,t_] := PDF[NormalDistribution[sigma1^2*t,sigma1*Sqrt[t]],x]

qy[x_,sigma1_,t_] := PDF[NormalDistribution[0,sigma1*Sqrt[t]],x]

qm[x_,mu_,sigma2_,t_] := PDF[NormalDistribution[mu*t,sigma2*Sqrt[t]],x]

pmax1[x_,k1_,sigma1_,t_,x0_,y0_]:= Max[qx[x,sigma1,t]*x0, k1*qy[x,sigma1,t]*y0]/(x0*(1-CDF[NormalDistribution[sigma1^2*t,sigma1*Sqrt[t]], sigma1^2*t/2-Log[x0/(k1*y0)]])+k1*y0(CDF[NormalDistribution[0,sigma1*Sqrt[t]], sigma1^2*t/2-Log[x0/(k1*y0)]]))

pmax1[x,k1,sigma1,t,x0,y0]

MATRIX

(1,1)

Simplify[Integrate[qx[x,sigma1,t]*qx[x,sigma1,t]/qm[x,mu,sigma2,t],{x,-Infinity,Infinity}],{sigma1>0,sigma2>0,t>=0}]

Simplify[Integrate[(qx[x,sigma1,t]/qm[x,mu,sigma2,t]-1)*(qx[x,sigma1,t]/qm[x,mu,sigma2,t]-1)*qm[x,mu,sigma2,t],{x,-Infinity,Infinity}],{sigma1>0,sigma2>0,t>=0}]

q11[mu_,t_,sigma1_, sigma2_] := (((E^(-(((((mu-(sigma1^2)))^2) t)/((sigma1^2)-2 (sigma2^2))))) sigma2)/((sigma1^2) (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))])))-1

(1,2)

Simplify[Integrate[qx[x,sigma1,t]*qy[x,sigma1,t]/qm[x,mu,sigma2,t],{x,-Infinity,Infinity}],{sigma1>0,sigma2>0,t>=0,sigma1<sigma2}]

q12[mu_,t_,sigma1_, sigma2_]:= (((E^(-(((2 (mu^2)-2 mu (sigma1^2)+(sigma1^4)-(sigma1^2) (sigma2^2)) t)/(2 ((sigma1^2)-2 (sigma2^2)))))) sigma2)/((sigma1^2) (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))])))-1

(2,2)

Simplify[Integrate[qy[x,sigma1,t]*qy[x,sigma1,t]/qm[x,mu,sigma2,t],{x,-Infinity,Infinity}],{sigma1>0,sigma2>0,t>=0,sigma1<sigma2}]

q22[mu_,t_,sigma1_, sigma2_]:=(((E^(-(((mu^2) t)/((sigma1^2)-2 (sigma2^2))))) (sigma2^2))/(sigma1 (Sqrt[-(sigma1^2)+2 (sigma2^2)])))-1

(1,3),(1,4)

Simplify[Integrate[(((((E^(-((x^2)/(2 (sigma1^2) t)))) k1 y0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))*qx[x,sigma1,t]/qm[x,mu,sigma2,t])/((1/2) k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))),{x,-Infinity,(t*sigma1^2)/2-Log[x0/(y0*k1)]}]+Integrate[(((((E^(-((((-(sigma1^2) t+x))^2)/(2 (sigma1^2) t)))) x0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))*qx[x,sigma1,t]/qm[x,mu,sigma2,t])/((1/2) k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))),{x,(t*sigma1^2)/2-Log[x0/(y0*k1)],Infinity}],{sigma1>0,t>0,k1>0,x0>0,y0>0}]

q13[mu_,k1_,sigma1_,sigma2_,t_,x0_,y0_]:=((E^(-(((2 (mu^2)-2 mu (sigma1^2)+(sigma1^4)-(sigma1^2) (sigma2^2)) t)/(2 ((sigma1^2)-2 (sigma2^2)))))) k1 sigma2 y0 ((1/(Sqrt[(((((mu-(sigma2^2)))^2) t)/(-(sigma1^2) (sigma2^2)+2 (sigma2^4)))]))sigma1 t ((mu-(sigma2^2)) Erf[sigma1 (Sqrt[(((((mu-(sigma2^2)))^2) t)/(-2 (sigma1^2) (sigma2^2)+4 (sigma2^4)))])]-sigma1 (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (sigma2^2) (Sqrt[((((mu-(sigma2^2)))^2)/(-(sigma1^2) (sigma2^2)+2 (sigma2^4)))]) (1+Erf[(((sigma1^2) (-mu+(sigma2^2)) (Sqrt[(t/(sigma1^2))-(t/(2 (sigma2^2)))]))/((sigma1^2)-2 (sigma2^2)))]))+((sigma1 Erf[((Sqrt[-(((((sigma1^2) (-2 mu+(sigma1^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))^2)/((sigma1^2) (sigma2^2) ((sigma1^2)-2 (sigma2^2)) t))])/(2 (Sqrt[2])))] ((sigma1^2) (-2 mu+(sigma1^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))/(Sqrt[(((((sigma1^2) (-2 mu+(sigma1^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))^2)/((-(sigma1^2) (sigma2^2)+2 (sigma2^4)) t))]))))/((sigma1^2) ((sigma1^2)-2 (sigma2^2)) (Sqrt[t]) (k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]-x0 (-2+Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])))-((E^(-(((((mu-(sigma1^2)))^2) t)/((sigma1^2)-2 (sigma2^2))))) sigma2 x0 (1+(((mu-2 (sigma2^2)) Erf[sigma1 (Sqrt[(((((mu-2 (sigma2^2)))^2) t)/(-2 (sigma1^2) (sigma2^2)+4 (sigma2^4)))])])/(sigma1 (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (sigma2^2) (Sqrt[((((mu-2 (sigma2^2)))^2)/(-(sigma1^2) (sigma2^2)+2 (sigma2^4)))])))-Erf[(((sigma1^2) (-mu+2 (sigma2^2)) (Sqrt[(t/(sigma1^2))-(t/(2 (sigma2^2)))]))/((sigma1^2)-2 (sigma2^2)))]+(Erf[((Sqrt[-(((((sigma1^2) (-2 mu+(sigma1^2)+2 (sigma2^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))^2)/((sigma1^2) (sigma2^2) ((sigma1^2)-2 (sigma2^2)) t))])/(2 (Sqrt[2])))] ((sigma1^2) (-2 mu+(sigma1^2)+2 (sigma2^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))/(sigma1 (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (sigma2^2) (Sqrt[(((((sigma1^2) (-2 mu+(sigma1^2)+2 (sigma2^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))^2)/(-(sigma1^2) (sigma2^2)+2 (sigma2^4)))]))))/((sigma1^2) (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (-k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (-2+Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])))-1

(2,3),(2,4)

Simplify[Integrate[(((((E^(-((x^2)/(2 (sigma1^2) t)))) k1 y0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))*qy[x,sigma1,t]/qm[x,mu,sigma2,t])/((1/2) k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))),{x,-Infinity,(t*sigma1^2)/2-Log[x0/(y0*k1)]}]+Integrate[(((((E^(-((((-(sigma1^2) t+x))^2)/(2 (sigma1^2) t)))) x0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))*qy[x,sigma1,t]/qm[x,mu,sigma2,t])/((1/2) k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))),{x,(t*sigma1^2)/2-Log[x0/(y0*k1)],Infinity}],{sigma1>0,t>0,k1>0,x0>0,y0>0}]

q23[mu_,k1_,sigma1_,sigma2_,t_,x0_,y0_]:=-(((E^(-(((2 (mu^2)-2 mu (sigma1^2)+(sigma1^4)-(sigma1^2) (sigma2^2)) t)/(2 ((sigma1^2)-2 (sigma2^2)))))) sigma2 x0 ((1/(Sqrt[(((((mu-(sigma2^2)))^2) t)/(-(sigma1^2) (sigma2^2)+2 (sigma2^4)))]))sigma1 t ((-mu+(sigma2^2)) Erf[sigma1 (Sqrt[(((((mu-(sigma2^2)))^2) t)/(-2 (sigma1^2) (sigma2^2)+4 (sigma2^4)))])]+sigma1 (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (sigma2^2) (Sqrt[((((mu-(sigma2^2)))^2)/(-(sigma1^2) (sigma2^2)+2 (sigma2^4)))]) (-1+Erf[(((sigma1^2) (-mu+(sigma2^2)) (Sqrt[(t/(sigma1^2))-(t/(2 (sigma2^2)))]))/((sigma1^2)-2 (sigma2^2)))]))+((sigma1 Erf[((Sqrt[-(((((sigma1^2) (-2 mu+(sigma1^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))^2)/((sigma1^2) (sigma2^2) ((sigma1^2)-2 (sigma2^2)) t))])/(2 (Sqrt[2])))] (-(sigma1^2) (-2 mu+(sigma1^2)) t+2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))/(Sqrt[(((((sigma1^2) (-2 mu+(sigma1^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))^2)/((-(sigma1^2) (sigma2^2)+2 (sigma2^4)) t))]))))/((sigma1^2) ((sigma1^2)-2 (sigma2^2)) (Sqrt[t]) (-k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (-2+Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))))+((E^(-(((mu^2) t)/((sigma1^2)-2 (sigma2^2))))) k1 sigma2 y0 (1-((mu Erf[sigma1 (Sqrt[(((mu^2) t)/(-2 (sigma1^2) (sigma2^2)+4 (sigma2^4)))])])/(sigma1 (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (sigma2^2) (Sqrt[((mu^2)/(-(sigma1^2) (sigma2^2)+2 (sigma2^4)))])))-Erf[((mu (sigma1^2) (Sqrt[(t/(sigma1^2))-(t/(2 (sigma2^2)))]))/((sigma1^2)-2 (sigma2^2)))]+(Erf[((Sqrt[-(((((sigma1^2) (2 mu-(sigma1^2)+2 (sigma2^2)) t+2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))^2)/((sigma1^2) (sigma2^2) ((sigma1^2)-2 (sigma2^2)) t))])/(2 (Sqrt[2])))] ((sigma1^2) (2 mu-(sigma1^2)+2 (sigma2^2)) t+2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))/(sigma1 (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (sigma2^2) (Sqrt[(((((sigma1^2) (2 mu-(sigma1^2)+2 (sigma2^2)) t+2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))^2)/(-(sigma1^2) (sigma2^2)+2 (sigma2^4)))]))))/((sigma1^2) (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]-x0 (-2+Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])))-1

(3,3),(4,4)

Simplify[Integrate[(((((E^(-((x^2)/(2 (sigma1^2) t)))) k1 y0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))/((1/2) k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))))^2/qm[x,mu,sigma2,t],{x,-Infinity,(t*sigma1^2)/2-Log[x0/(y0*k1)]}]+Integrate[(((((E^(-((((-(sigma1^2) t+x))^2)/(2 (sigma1^2) t)))) x0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))/((1/2) k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))))^2/qm[x,mu,sigma2,t],{x,(t*sigma1^2)/2-Log[x0/(y0*k1)],Infinity}],{sigma1>0,t>0,k1>0,x0>0,y0>0}]

q33[mu_,k1_,sigma1_,sigma2_,t_,x0_,y0_]:=(2 sigma2 ((E^(-(((((mu-(sigma1^2)))^2) t)/((sigma1^2)-2 (sigma2^2))))) (x0^2) (1+(((mu-2 (sigma2^2)) Erf[sigma1 (Sqrt[(((((mu-2 (sigma2^2)))^2) t)/(-2 (sigma1^2) (sigma2^2)+4 (sigma2^4)))])])/(sigma1 (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (sigma2^2) (Sqrt[((((mu-2 (sigma2^2)))^2)/(-(sigma1^2) (sigma2^2)+2 (sigma2^4)))])))-Erf[(((sigma1^2) (-mu+2 (sigma2^2)) (Sqrt[(t/(sigma1^2))-(t/(2 (sigma2^2)))]))/((sigma1^2)-2 (sigma2^2)))]+(Erf[((Sqrt[-(((((sigma1^2) (-2 mu+(sigma1^2)+2 (sigma2^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))^2)/((sigma1^2) (sigma2^2) ((sigma1^2)-2 (sigma2^2)) t))])/(2 (Sqrt[2])))] ((sigma1^2) (-2 mu+(sigma1^2)+2 (sigma2^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))/(sigma1 (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (sigma2^2) (Sqrt[(((((sigma1^2) (-2 mu+(sigma1^2)+2 (sigma2^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))^2)/(-(sigma1^2) (sigma2^2)+2 (sigma2^4)))])))+(E^(-(((mu^2) t)/((sigma1^2)-2 (sigma2^2))))) (k1^2) (y0^2) (1-((mu Erf[sigma1 (Sqrt[(((mu^2) t)/(-2 (sigma1^2) (sigma2^2)+4 (sigma2^4)))])])/(sigma1 (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (sigma2^2) (Sqrt[((mu^2)/(-(sigma1^2) (sigma2^2)+2 (sigma2^4)))])))-Erf[((mu (sigma1^2) (Sqrt[(t/(sigma1^2))-(t/(2 (sigma2^2)))]))/((sigma1^2)-2 (sigma2^2)))]+(Erf[((Sqrt[-(((((sigma1^2) (2 mu-(sigma1^2)+2 (sigma2^2)) t+2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))^2)/((sigma1^2) (sigma2^2) ((sigma1^2)-2 (sigma2^2)) t))])/(2 (Sqrt[2])))] ((sigma1^2) (2 mu-(sigma1^2)+2 (sigma2^2)) t+2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))/(sigma1 (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (sigma2^2) (Sqrt[(((((sigma1^2) (2 mu-(sigma1^2)+2 (sigma2^2)) t+2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))^2)/(-(sigma1^2) (sigma2^2)+2 (sigma2^4)))])))))/((sigma1^2) (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (((k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]-x0 (-2+Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])))^2))-1

(3,4) WLOG and for the sake of the calculation we will assume .

Simplify[Integrate[((((E^(-((x^2)/(2 (sigma1^2) t)))) k1 y0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))/((1/2) k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])))*((((E^(-((x^2)/(2 (sigma1^2) t)))) k2 y0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))/((1/2) k2 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])))/qm[x,mu,sigma2,t],{x,-Infinity,(t*sigma1^2)/2-Log[x0/(y0*k1)]},Assumptions-> k1<k2]+Integrate[(((((E^(-((((-(sigma1^2) t+x))^2)/(2 (sigma1^2) t)))) x0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))/((1/2) k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))))*((((E^(-((x^2)/(2 (sigma1^2) t)))) k2 y0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))/((1/2) k2 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])))/qm[x,mu,sigma2,t],{x,(t*sigma1^2)/2-Log[x0/(y0*k1)],(t*sigma1^2)/2-Log[x0/(y0*k2)]},Assumptions-> k1<k2]+Integrate[(((((E^(-((((-(sigma1^2) t+x))^2)/(2 (sigma1^2) t)))) x0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))/((1/2) k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))))*(((((E^(-((((-(sigma1^2) t+x))^2)/(2 (sigma1^2) t)))) x0)/((Sqrt[2 π]) sigma1 (Sqrt[t])))/((1/2) k2 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]+x0 (1-(1/2) Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]))))/qm[x,mu,sigma2,t],{x,(t*sigma1^2)/2-Log[x0/(y0*k2)],Infinity},Assumptions-> k1<k2],{sigma1>0,t>0,k1>0,x0>0,y0>0, k1<k2}]

q34[mu_,k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_]:=(2 sigma2 ((1/(Sqrt[-((-(2/(sigma1^2))+(1/(sigma2^2)))/t)]))(E^(-(((mu^2) t)/((sigma1^2)-2 (sigma2^2))))) k1 k2 (y0^2) (1-((mu Erf[sigma1 (Sqrt[(((mu^2) t)/(-2 (sigma1^2) (sigma2^2)+4 (sigma2^4)))])])/(sigma1 (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (sigma2^2) (Sqrt[((mu^2)/(-(sigma1^2) (sigma2^2)+2 (sigma2^4)))])))-Erf[((mu (sigma1^2) (Sqrt[(t/(sigma1^2))-(t/(2 (sigma2^2)))]))/((sigma1^2)-2 (sigma2^2)))]+(Erf[((Sqrt[-(((((sigma1^2) (2 mu-(sigma1^2)+2 (sigma2^2)) t+2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))^2)/((sigma1^2) (sigma2^2) ((sigma1^2)-2 (sigma2^2)) t))])/(2 (Sqrt[2])))] ((sigma1^2) (2 mu-(sigma1^2)+2 (sigma2^2)) t+2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))/(sigma1 (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (sigma2^2) (Sqrt[(((((sigma1^2) (2 mu-(sigma1^2)+2 (sigma2^2)) t+2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))^2)/(-(sigma1^2) (sigma2^2)+2 (sigma2^4)))])))+(1/(Sqrt[-((-(2/(sigma1^2))+(1/(sigma2^2)))/t)]))(E^(-(((((mu-(sigma1^2)))^2) t)/((sigma1^2)-2 (sigma2^2))))) (x0^2) (1+(((mu-2 (sigma2^2)) Erf[sigma1 (Sqrt[(((((mu-2 (sigma2^2)))^2) t)/(-2 (sigma1^2) (sigma2^2)+4 (sigma2^4)))])])/(sigma1 (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (sigma2^2) (Sqrt[((((mu-2 (sigma2^2)))^2)/(-(sigma1^2) (sigma2^2)+2 (sigma2^4)))])))-Erf[(((sigma1^2) (-mu+2 (sigma2^2)) (Sqrt[(t/(sigma1^2))-(t/(2 (sigma2^2)))]))/((sigma1^2)-2 (sigma2^2)))]+(Erf[((Sqrt[-(((((sigma1^2) (-2 mu+(sigma1^2)+2 (sigma2^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k2 y0))]))^2)/((sigma1^2) (sigma2^2) ((sigma1^2)-2 (sigma2^2)) t))])/(2 (Sqrt[2])))] ((sigma1^2) (-2 mu+(sigma1^2)+2 (sigma2^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k2 y0))]))/(sigma1 (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) (sigma2^2) (Sqrt[(((((sigma1^2) (-2 mu+(sigma1^2)+2 (sigma2^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k2 y0))]))^2)/(-(sigma1^2) (sigma2^2)+2 (sigma2^4)))])))+(1/(-(sigma1^2)+2 (sigma2^2)))(E^(-(((2 (mu^2)-2 mu (sigma1^2)+(sigma1^4)-(sigma1^2) (sigma2^2)) t)/(2 ((sigma1^2)-2 (sigma2^2)))))) k2 x0 y0 (((sigma1 Erf[((Sqrt[-(((((sigma1^2) (-2 mu+(sigma1^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))^2)/((sigma1^2) (sigma2^2) ((sigma1^2)-2 (sigma2^2)) t))])/(2 (Sqrt[2])))] ((sigma1^2) (-2 mu+(sigma1^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))/(Sqrt[(((((sigma1^2) (-2 mu+(sigma1^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k1 y0))]))^2)/((-(sigma1^2) (sigma2^2)+2 (sigma2^4)) t))]))+((sigma1 Erf[((Sqrt[-(((((sigma1^2) (-2 mu+(sigma1^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k2 y0))]))^2)/((sigma1^2) (sigma2^2) ((sigma1^2)-2 (sigma2^2)) t))])/(2 (Sqrt[2])))] (-(sigma1^2) (-2 mu+(sigma1^2)) t+2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k2 y0))]))/(Sqrt[(((((sigma1^2) (-2 mu+(sigma1^2)) t-2 ((sigma1^2)-2 (sigma2^2)) Log[(x0/(k2 y0))]))^2)/((-(sigma1^2) (sigma2^2)+2 (sigma2^4)) t))])))))/((sigma1^2) (Sqrt[t]) (k1 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]-x0 (-2+Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k1 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])) (k2 y0 Erfc[((-(((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))]-x0 (-2+Erfc[(((((sigma1^2) t)/2)+Log[(x0/(k2 y0))])/((Sqrt[2]) sigma1 (Sqrt[t])))])))-1

No contract portfolio

Find the weights:

A0[sigma1_,sigma2_,t_,x0_,y0_,mu_]:= {{q11[mu,t,sigma1, sigma2],q12[mu,t,sigma1, sigma2]},{q12[mu,t,sigma1, sigma2],q22[mu,t,sigma1, sigma2]}}

A0[0.3,0.3,10,1,1,0.05]

A:={{q11,q12},{q12,q22}}

Solve[D[{w1,1-w1}.A0[sigma1,sigma2,t,x0,y0,mu].{{w1},{1-w1}} ,w1]==0 
, {w1}]

Simplify[((E^(((((mu-(sigma1^2)))^2) t)/((sigma1^2)-2 (sigma2^2)))) ((E^(((2 (mu^2)-2 mu (sigma1^2)+(sigma1^4)-(sigma1^2) (sigma2^2)) t)/(2 ((sigma1^2)-2 (sigma2^2))))) sigma1 (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) sigma2-(E^(((mu^2) t)/((sigma1^2)-2 (sigma2^2)))) (Sqrt[-(sigma1^2)+2 (sigma2^2)])))/((E^((((((mu-(sigma1^2)))^2) t)/((sigma1^2)-2 (sigma2^2)))+(((2 (mu^2)-2 mu (sigma1^2)+(sigma1^4)-(sigma1^2) (sigma2^2)) t)/(2 ((sigma1^2)-2 (sigma2^2)))))) sigma1 (Sqrt[(2/(sigma1^2))-(1/(sigma2^2))]) sigma2-2 (E^((((mu^2) t)/((sigma1^2)-2 (sigma2^2)))+(((((mu-(sigma1^2)))^2) t)/((sigma1^2)-2 (sigma2^2))))) (Sqrt[-(sigma1^2)+2 (sigma2^2)])+(E^((((mu^2) t)/((sigma1^2)-2 (sigma2^2)))+(((2 (mu^2)-2 mu (sigma1^2)+(sigma1^4)-(sigma1^2) (sigma2^2)) t)/(2 ((sigma1^2)-2 (sigma2^2)))))) (Sqrt[-(sigma1^2)+2 (sigma2^2)])),{sigma1>0,sigma2>0}]

Tw[sigma1_,sigma2_,t_,mu_]:= (((E^(((sigma1^4) t)/((sigma1^2)-2 (sigma2^2))))-(E^(((sigma1^2) (2 mu+(sigma1^2)+(sigma2^2)) t)/(2 ((sigma1^2)-2 (sigma2^2))))))/((E^((2 mu (sigma1^2) t)/((sigma1^2)-2 (sigma2^2))))+(E^(((sigma1^4) t)/((sigma1^2)-2 (sigma2^2))))-2 (E^(((sigma1^2) (2 mu+(sigma1^2)+(sigma2^2)) t)/(2 ((sigma1^2)-2 (sigma2^2)))))))

Tw[0.2,0.2,10,0.03]

Calculate the payoff:

TunedPayoff0[sigma1_,sigma2_,t_,mu_]:= NIntegrate[qm[x,mu,sigma2,t]*Log[(qx[x,sigma1,t]*Tw[sigma1,sigma2,t,mu]+(1-Tw[sigma1,sigma2,t,mu])*qy[x,sigma1,t])/qy[x,sigma1,t]],{x,-Infinity,Infinity}]

TunedDensity0[x_,mu_,sigma1_,sigma2_,t_]:=(qx[x,sigma1,t]*Tw[sigma1,sigma2,t,mu]+(1-Tw[sigma1,sigma2,t,mu])*qy[x,sigma1,t])

Tw[0.3,0.3,10,0.05]

TunedPayoff0[0.2,0.2,10,0.03]

1 contract portfolio

Find the weights:

A:= {{q11,q12,q13},{q12,q22,q23},{q13,q23,q33}}

Solve[D[{w1,w2,1-w1-w2}.A.{{w1},{w2},{1-w1-w2}},w1]==0&& D[{w1,w2,1-w1-w2}.A.{{w1},{w2},{1-w1-w2}},w2]==0, {w1,w2}]

Tw1[mu_,k1_,sigma1_,sigma2_,t_,x0_,y0_]:= -((-q13[mu,k1,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2]+q12[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0]+q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0]-((q23[mu,k1,sigma1,sigma2,t,x0,y0])^2)-q12[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0]+q22[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0])/(((q12[mu,t,sigma1, sigma2])^2)-2 q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0]+((q13[mu,k1,sigma1,sigma2,t,x0,y0])^2)-q11[mu,t,sigma1, sigma2] q22[mu,t,sigma1, sigma2]+2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2]+2 q11[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0]-2 q12[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0]-2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0]+((q23[mu,k1,sigma1,sigma2,t,x0,y0])^2)-q11[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0]-q22[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0]))

Tw1[0.03,2,0.2,0.2,20,1,1]

Tw2[mu_,k1_,sigma1_,sigma2_,t_,x0_,y0_]:= -((q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0]-((q13[mu,k1,sigma1,sigma2,t,x0,y0])^2)-q11[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0]+q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0]+q11[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0]-q12[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0])/(((q12[mu,t,sigma1, sigma2])^2)-2 q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0]+((q13[mu,k1,sigma1,sigma2,t,x0,y0])^2)-q11[mu,t,sigma1, sigma2] q22[mu,t,sigma1, sigma2]+2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2]+2 q11[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0]-2 q12[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0]-2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0]+((q23[mu,k1,sigma1,sigma2,t,x0,y0])^2)-q11[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0]-q22[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0]))

Calculate payoff and optimal strike:

TunedPayoff1[mu_,sigma1_,sigma2_,t_,k_,x0_,y0_]:= NIntegrate[Log[(Tw1[mu,k,sigma1,sigma2,t,x0,y0]*qx[x,sigma1,t] + Tw2[mu,k,sigma1,sigma2,t,x0,y0]qy[x,sigma1,t] + (1-Tw1[mu,k,sigma1,sigma2,t,x0,y0]-Tw2[mu,k,sigma1,sigma2,t,x0,y0])*pmax1[x,k,sigma1,t,x0,y0])/qy[x,sigma1,t]]*qm[x,mu,sigma2,t],{x,-Infinity,Infinity}]

TunedDensity1[x_,mu_,sigma1_,sigma2_,t_,k_,x0_,y0_]:= Tw1[mu,k,sigma1,sigma2,t,x0,y0]*qx[x,sigma1,t] + Tw2[mu,k,sigma1,sigma2,t,x0,y0]qy[x,sigma1,t] + (1-Tw1[mu,k,sigma1,sigma2,t,x0,y0]-Tw2[mu,k,sigma1,sigma2,t,x0,y0])*pmax1[x,k,sigma1,t,x0,y0]

Tw2[0.03,2,0.2,0.2,20,1,1]

BT1[mu_,k1_,sigma1_,sigma2_,t_,x0_,y0_]:= {{q11[mu,t,sigma1, sigma2],q12[mu,t,sigma1, sigma2],q13[mu,k1,sigma1,sigma2,t,x0,y0]},{q12[mu,t,sigma1, sigma2],q22[mu,t,sigma1, sigma2],q23[mu,k1,sigma1,sigma2,t,x0,y0]},{q13[mu,k1,sigma1,sigma2,t,x0,y0],q23[mu,k1,sigma1,sigma2,t,x0,y0],q33[mu,k1,sigma1,sigma2,t,x0,y0]}}

TunedApMax1[mu_,sigma1_,sigma2_,t_,x0_,y0_]:= FindMinimum[{Norm[{Tw1[mu,k,sigma1,sigma2,t,x0,y0],Tw2[mu,k,sigma1,sigma2,t,x0,y0],1-Tw1[mu,k,sigma1,sigma2,t,x0,y0]-Tw2[mu,k,sigma1,sigma2,t,x0,y0]}.BT1[mu,k,sigma1,sigma2,t,x0,y0].{{Tw1[mu,k,sigma1,sigma2,t,x0,y0]},{Tw2[mu,k,sigma1,sigma2,t,x0,y0]},{1-Tw1[mu,k,sigma1,sigma2,t,x0,y0]-Tw2[mu,k,sigma1,sigma2,t,x0,y0]}}], k>0},{k}]

TunedApMaxPlot1[mu_,sigma1_,sigma2_,t_,x0_,y0_]:= Plot[{Norm[{Tw1[mu,k,sigma1,sigma2,t,x0,y0],Tw2[mu,k,sigma1,sigma2,t,x0,y0],1-Tw1[mu,k,sigma1,sigma2,t,x0,y0]-Tw2[mu,k,sigma1,sigma2,t,x0,y0]}.BT1[mu,k,sigma1,sigma2,t,x0,y0].{{Tw1[mu,k,sigma1,sigma2,t,x0,y0]},{Tw2[mu,k,sigma1,sigma2,t,x0,y0]},{1-Tw1[mu,k,sigma1,sigma2,t,x0,y0]-Tw2[mu,k,sigma1,sigma2,t,x0,y0]}}]},{k,0,5}]

TunedApMax1[0.05,0.3,0.3,10,1,1]

TunedApMax1[0.05,0.3,0.22,10,1,1]

TunedApMax1[0.05,0.3,0.35,1,1,1]

Tw1[0.05,1.173,0.3,0.22,10,1,1]

Tw2[0.05,1.173,0.3,0.22,10,1,1]

Tw2[0.05,1.0193,0.3,0.3,10,1,1]

Tw1[0.05,1.0193,0.3,0.3,10,1,1]

1-Tw1[0.05,1.0193,0.3,0.3,10,1,1]-Tw2[0.05,1.0193,0.3,0.3,10,1,1]

TunedPayoff1[0.05,0.3,0.3,10,1.0193,1,1]

One contract option strategy:

STRADDLE: we bet on larger volatility than expected by the market, i. e. .

Define option prices according to Black-Scholes model:

CallOptionPrice[sigma1_,t_,xy_,k_]:= xy*CDF[NormalDistribution[0,1], sigma1*Sqrt[t]/2+1/(sigma1*Sqrt[t])*Log[xy/k]] - k*CDF[NormalDistribution[0,1], -sigma1*Sqrt[t]/2+1/(sigma1*Sqrt[t])*Log[xy/k]]

PutOptionPrice[sigma1_,t_,xy_,k_]:= CallOptionPrice[sigma1,t,xy,k]+k-xy

StraddlePayoff[p_,sigma1_,t_,xy_,k_]:= -CallOptionPrice[sigma1,t,xy,k]*Boole[p<k] + Boole[p>=k]*(p-CallOptionPrice[sigma1,t,xy,k]-k) + Boole[p<k]*(-p+k-PutOptionPrice[sigma1,t,xy,k]) -PutOptionPrice[sigma1,t,xy,k]* Boole[p>=k]

CallOptionPrice[0.3,10,1,1.173]

PutOptionPrice[0.3,10,1,1.173]

TunedApMax1[0.05,0.3,0.22,10,1,1]

TunedApMax1[-0.05,0.3,0.3,10,1,1]

Plot[StraddlePayoff[p,0.3,10,1,1.173],{p,0,2.5},PlotLabel->"P/L of Straddle Strategy with optimal strike",AxesLabel->{"price","P/L"}]

2 contracts portfolio

Find the weights:

matice:= {{q11,q12,q13,q13k},{q12,q22,q23,q23k},{q13,q23,q33,q34},{q13k,q23k,q34,q33k}}

Solve[D[{w1,w2,w3,1-w1-w2-w3}.matice.{{w1},{w2},{w3},{1-w1-w2-w3}},w1]==0&& D[{w1,w2,w3,1-w1-w2-w3}.matice.{{w1},{w2},{w3},{1-w1-w2-w3}},w2]==0 &&
D[{w1,w2,w3,1-w1-w2-w3}.matice.{{w1},{w2},{w3},{1-w1-w2-w3}},w3]==0, {w1,w2,w3}]

v1[mu_,k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_]:= -((((2 q22[mu,t,sigma1, sigma2]-4 q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]) (2 q33[mu,k1,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]-4 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])-(((2 q23[mu,k1,sigma1,sigma2,t,x0,y0]-2 q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]))^2)) ((2 q13[mu,k2,sigma1,sigma2,t,x0,y0]-2 q33[mu,k2,sigma1,sigma2,t,x0,y0]) (2 q33[mu,k1,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]-4 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])-(2 q13[mu,k1,sigma1,sigma2,t,x0,y0]-2 q13[mu,k2,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]) (-2 q33[mu,k2,sigma1,sigma2,t,x0,y0]+2 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]))-((2 q12[mu,t,sigma1, sigma2]-2 q13[mu,k2,sigma1,sigma2,t,x0,y0]-2 q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]) (2 q33[mu,k1,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]-4 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])-(2 q13[mu,k1,sigma1,sigma2,t,x0,y0]-2 q13[mu,k2,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]) (2 q23[mu,k1,sigma1,sigma2,t,x0,y0]-2 q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])) ((2 q23[mu,k2,sigma1,sigma2,t,x0,y0]-2 q33[mu,k2,sigma1,sigma2,t,x0,y0]) (2 q33[mu,k1,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]-4 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])-(2 q23[mu,k1,sigma1,sigma2,t,x0,y0]-2 q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]) (-2 q33[mu,k2,sigma1,sigma2,t,x0,y0]+2 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])))/(-((((2 q12[mu,t,sigma1, sigma2]-2 q13[mu,k2,sigma1,sigma2,t,x0,y0]-2 q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]) (2 q33[mu,k1,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]-4 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])-(2 q13[mu,k1,sigma1,sigma2,t,x0,y0]-2 q13[mu,k2,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]) (2 q23[mu,k1,sigma1,sigma2,t,x0,y0]-2 q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])))^2)+((2 q11[mu,t,sigma1, sigma2]-4 q13[mu,k2,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]) (2 q33[mu,k1,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]-4 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])-(((2 q13[mu,k1,sigma1,sigma2,t,x0,y0]-2 q13[mu,k2,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]))^2)) ((2 q22[mu,t,sigma1, sigma2]-4 q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]) (2 q33[mu,k1,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]-4 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])-(((2 q23[mu,k1,sigma1,sigma2,t,x0,y0]-2 q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]))^2))))

v2[mu_,k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_]:= -((-q13[mu,k1,sigma1,sigma2,t,x0,y0] q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0]+((q13[mu,k2,sigma1,sigma2,t,x0,y0])^2) q23[mu,k1,sigma1,sigma2,t,x0,y0]+((q13[mu,k1,sigma1,sigma2,t,x0,y0])^2) q23[mu,k2,sigma1,sigma2,t,x0,y0]-q13[mu,k1,sigma1,sigma2,t,x0,y0] q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]+q12[mu,t,sigma1, sigma2] q13[mu,k2,sigma1,sigma2,t,x0,y0] q33[mu,k1,sigma1,sigma2,t,x0,y0]-((q13[mu,k2,sigma1,sigma2,t,x0,y0])^2) q33[mu,k1,sigma1,sigma2,t,x0,y0]-q11[mu,t,sigma1, sigma2] q23[mu,k2,sigma1,sigma2,t,x0,y0] q33[mu,k1,sigma1,sigma2,t,x0,y0]+q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0] q33[mu,k1,sigma1,sigma2,t,x0,y0]+q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]-((q13[mu,k1,sigma1,sigma2,t,x0,y0])^2) q33[mu,k2,sigma1,sigma2,t,x0,y0]-q11[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]+q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]+q11[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]-q12[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]-q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-q12[mu,t,sigma1, sigma2] q13[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q13[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+q11[mu,t,sigma1, sigma2]q23[mu,k1,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+q11[mu,t,sigma1, sigma2] q23[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-q11[mu,t,sigma1, sigma2]*((q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2)+q12[mu,t,sigma1, sigma2]*((q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2))/(((q13[mu,k1,sigma1,sigma2,t,x0,y0])^2) q22[mu,t,sigma1, sigma2]-2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q13[mu,k2,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2]+((q13[mu,k2,sigma1,sigma2,t,x0,y0])^2) q22[mu,t,sigma1, sigma2]-2 q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0]+2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0]-2 ((q13[mu,k2,sigma1,sigma2,t,x0,y0])^2) q23[mu,k1,sigma1,sigma2,t,x0,y0]+q11[mu,t,sigma1, sigma2] ((q23[mu,k1,sigma1,sigma2,t,x0,y0])^2)-2 q13[mu,k2,sigma1,sigma2,t,x0,y0] ((q23[mu,k1,sigma1,sigma2,t,x0,y0])^2)+2 q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]-2 ((q13[mu,k1,sigma1,sigma2,t,x0,y0])^2) q23[mu,k2,sigma1,sigma2,t,x0,y0]-2 q12[mu,t,sigma1, sigma2] q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q13[mu,k1,sigma1,sigma2,t,x0,y0]q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]-2 q11[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0]q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]+q11[mu,t,sigma1, sigma2] ((q23[mu,k2,sigma1,sigma2,t,x0,y0])^2)-2 q13[mu,k1,sigma1,sigma2,t,x0,y0] ((q23[mu,k2,sigma1,sigma2,t,x0,y0])^2)+((q12[mu,t,sigma1, sigma2])^2) q33[mu,k1,sigma1,sigma2,t,x0,y0]-2 q12[mu,t,sigma1, sigma2] q13[mu,k2,sigma1,sigma2,t,x0,y0] q33[mu,k1,sigma1,sigma2,t,x0,y0]+((q13[mu,k2,sigma1,sigma2,t,x0,y0])^2) q33[mu,k1,sigma1,sigma2,t,x0,y0]-q11[mu,t,sigma1, sigma2] q22[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0]+2 q13[mu,k2,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0]+2 q11[mu,t,sigma1, sigma2] q23[mu,k2,sigma1,sigma2,t,x0,y0] q33[mu,k1,sigma1,sigma2,t,x0,y0]-2 q12[mu,t,sigma1, sigma2] q23[mu,k2,sigma1,sigma2,t,x0,y0] q33[mu,k1,sigma1,sigma2,t,x0,y0]-2 q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0] q33[mu,k1,sigma1,sigma2,t,x0,y0]+((q23[mu,k2,sigma1,sigma2,t,x0,y0])^2) q33[mu,k1,sigma1,sigma2,t,x0,y0]+((q12[mu,t,sigma1, sigma2])^2) q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]+((q13[mu,k1,sigma1,sigma2,t,x0,y0])^2) q33[mu,k2,sigma1,sigma2,t,x0,y0]-q11[mu,t,sigma1, sigma2] q22[mu,t,sigma1, sigma2] q33[mu,k2,sigma1,sigma2,t,x0,y0]+2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2] q33[mu,k2,sigma1,sigma2,t,x0,y0]+2 q11[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 q12[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]+((q23[mu,k1,sigma1,sigma2,t,x0,y0])^2) q33[mu,k2,sigma1,sigma2,t,x0,y0]-q11[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]-q22[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 ((q12[mu,t,sigma1, sigma2])^2) q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q13[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q13[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 q11[mu,t,sigma1, sigma2] q22[mu,t,sigma1, sigma2] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 q13[mu,k2,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 q11[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 q11[mu,t,sigma1, sigma2] q23[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q23[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 q23[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+q11[mu,t,sigma1, sigma2] ((q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2)-2 q12[mu,t,sigma1, sigma2] ((q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2)+q22[mu,t,sigma1, sigma2] ((q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2)))

v3[mu_,k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_]:=-((q13[mu,k1,sigma1,sigma2,t,x0,y0] q13[mu,k2,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2]-((q13[mu,k2,sigma1,sigma2,t,x0,y0])^2) q22[mu,t,sigma1, sigma2]-q12[mu,t,sigma1, sigma2] q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0]+((q13[mu,k2,sigma1,sigma2,t,x0,y0])^2) q23[mu,k1,sigma1,sigma2,t,x0,y0]-q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]-q13[mu,k1,sigma1,sigma2,t,x0,y0] q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]+q11[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]-q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]-q11[mu,t,sigma1, sigma2] ((q23[mu,k2,sigma1,sigma2,t,x0,y0])^2)+q13[mu,k1,sigma1,sigma2,t,x0,y0] ((q23[mu,k2,sigma1,sigma2,t,x0,y0])^2)-((q12[mu,t,sigma1, sigma2])^2) q33[mu,k2,sigma1,sigma2,t,x0,y0]+q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]+q11[mu,t,sigma1, sigma2] q22[mu,t,sigma1, sigma2] q33[mu,k2,sigma1,sigma2,t,x0,y0]-q13[mu,k1,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2] q33[mu,k2,sigma1,sigma2,t,x0,y0]-q11[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]+q12[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]+((q12[mu,t,sigma1, sigma2])^2) q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-q12[mu,t,sigma1, sigma2] q13[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-q11[mu,t,sigma1, sigma2] q22[mu,t,sigma1, sigma2] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+q13[mu,k2,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+q11[mu,t,sigma1, sigma2] q23[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-q12[mu,t,sigma1, sigma2] q23[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])/(((q13[mu,k1,sigma1,sigma2,t,x0,y0])^2) q22[mu,t,sigma1, sigma2]-2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q13[mu,k2,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2]+((q13[mu,k2,sigma1,sigma2,t,x0,y0])^2) q22[mu,t,sigma1, sigma2]-2 q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0]+2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0]-2 ((q13[mu,k2,sigma1,sigma2,t,x0,y0])^2) q23[mu,k1,sigma1,sigma2,t,x0,y0]+q11[mu,t,sigma1, sigma2] ((q23[mu,k1,sigma1,sigma2,t,x0,y0])^2)-2 q13[mu,k2,sigma1,sigma2,t,x0,y0] ((q23[mu,k1,sigma1,sigma2,t,x0,y0])^2)+2 q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]-2 ((q13[mu,k1,sigma1,sigma2,t,x0,y0])^2) q23[mu,k2,sigma1,sigma2,t,x0,y0]-2 q12[mu,t,sigma1, sigma2] q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]-2 q11[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]+q11[mu,t,sigma1, sigma2] ((q23[mu,k2,sigma1,sigma2,t,x0,y0])^2)-2 q13[mu,k1,sigma1,sigma2,t,x0,y0] ((q23[mu,k2,sigma1,sigma2,t,x0,y0])^2)+((q12[mu,t,sigma1, sigma2])^2) q33[mu,k1,sigma1,sigma2,t,x0,y0]-2 q12[mu,t,sigma1, sigma2] q13[mu,k2,sigma1,sigma2,t,x0,y0] q33[mu,k1,sigma1,sigma2,t,x0,y0]+((q13[mu,k2,sigma1,sigma2,t,x0,y0])^2) q33[mu,k1,sigma1,sigma2,t,x0,y0]-q11[mu,t,sigma1, sigma2] q22[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0]+2 q13[mu,k2,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0]+2 q11[mu,t,sigma1, sigma2] q23[mu,k2,sigma1,sigma2,t,x0,y0] q33[mu,k1,sigma1,sigma2,t,x0,y0]-2 q12[mu,t,sigma1, sigma2] q23[mu,k2,sigma1,sigma2,t,x0,y0] q33[mu,k1,sigma1,sigma2,t,x0,y0]-2 q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0] q33[mu,k1,sigma1,sigma2,t,x0,y0]+((q23[mu,k2,sigma1,sigma2,t,x0,y0])^2) q33[mu,k1,sigma1,sigma2,t,x0,y0]+((q12[mu,t,sigma1, sigma2])^2) q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]+((q13[mu,k1,sigma1,sigma2,t,x0,y0])^2) q33[mu,k2,sigma1,sigma2,t,x0,y0]-q11[mu,t,sigma1, sigma2] q22[mu,t,sigma1, sigma2] q33[mu,k2,sigma1,sigma2,t,x0,y0]+2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2] q33[mu,k2,sigma1,sigma2,t,x0,y0]+2 q11[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 q12[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]+((q23[mu,k1,sigma1,sigma2,t,x0,y0])^2) q33[mu,k2,sigma1,sigma2,t,x0,y0]-q11[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]-q22[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0] q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 ((q12[mu,t,sigma1, sigma2])^2) q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q13[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q13[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 q11[mu,t,sigma1, sigma2] q22[mu,t,sigma1, sigma2] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 q13[mu,k2,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 q11[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 q11[mu,t,sigma1, sigma2] q23[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q23[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 q23[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+q11[mu,t,sigma1, sigma2] ((q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2)-2 q12[mu,t,sigma1, sigma2] ((q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2)+q22[mu,t,sigma1, sigma2] ((q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])^2)))

Find the optimal strikes and payoffs

matka1[mu_,k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_]:= {{q11[mu,t,sigma1, sigma2],q12[mu,t,sigma1, sigma2],q13[mu,k1,sigma1,sigma2,t,x0,y0],q13[mu,k2,sigma1,sigma2,t,x0,y0]},{q12[mu,t,sigma1, sigma2],q22[mu,t,sigma1, sigma2],q23[mu,k1,sigma1,sigma2,t,x0,y0],q23[mu,k2,sigma1,sigma2,t,x0,y0]},{q13[mu,k1,sigma1,sigma2,t,x0,y0],q23[mu,k1,sigma1,sigma2,t,x0,y0],q33[mu,k1,sigma1,sigma2,t,x0,y0],q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]},{q13[mu,k2,sigma1,sigma2,t,x0,y0],q23[mu,k2,sigma1,sigma2,t,x0,y0],q34[mu,k1,k2,sigma1,sigma2,t,x0,y0],q33[mu,k2,sigma1,sigma2,t,x0,y0]}}

TunedApMax2[mu_,sigma1_,sigma2_,t_,x0_,y0_]:= FindMinimum[{Norm[{v1[mu,k1,k2,sigma1,sigma2,t,x0,y0],v2[mu,k1,k2,sigma1,sigma2,t,x0,y0],v3[mu,k1,k2,sigma1,sigma2,t,x0,y0],1-v1[mu,k1,k2,sigma1,sigma2,t,x0,y0]-v2[mu,k1,k2,sigma1,sigma2,t,x0,y0]-v3[mu,k1,k2,sigma1,sigma2,t,x0,y0]}.matka1[mu,k1,k2,sigma1,sigma2,t,x0,y0].{{v1[mu,k1,k2,sigma1,sigma2,t,x0,y0]},{v2[mu,k1,k2,sigma1,sigma2,t,x0,y0]},{v3[mu,k1,k2,sigma1,sigma2,t,x0,y0]},{1-v1[mu,k1,k2,sigma1,sigma2,t,x0,y0]-v2[mu,k1,k2,sigma1,sigma2,t,x0,y0]-v3[mu,k1,k2,sigma1,sigma2,t,x0,y0]}}],0.1<k1<k2<10},{k1,k2}]

TunedPayoff2[mu_,k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_]:= NIntegrate[Log[(v1[mu,k1,k2,sigma1,sigma2,t,x0,y0]*qx[x,sigma1,t] + v2[mu,k1,k2,sigma1,sigma2,t,x0,y0]qy[x,sigma1,t] + v3[mu,k1,k2,sigma1,sigma2,t,x0,y0]*pmax1[x,k1,sigma1,t,x0,y0]+ (1-v1[mu,k1,k2,sigma1,sigma2,t,x0,y0]-v2[mu,k1,k2,sigma1,sigma2,t,x0,y0]-v3[mu,k1,k2,sigma1,sigma2,t,x0,y0])*pmax1[x,k2,sigma1,t,x0,y0])/qy[x,sigma1,t]]*qm[x,mu,sigma2,t],{x,-Infinity,Infinity}]

TunedApMax2[0.05,0.3,0.4,10,1,1]

v1[0.05,0.27936,4.0873,0.3,0.4,10,1,1]

v2[0.05,0.27936,4.0873,0.3,0.4,10,1,1]

v3[0.05,0.27936,4.0873,0.3,0.4,10,1,1]

TunedPayoff2[0.05,0.27936,4.0873,0.3,0.4,10,1,1]

completeDKL[0.05,0.3,0.4,10]

OPTIONS HAVE THE SAME WEIGHTS:

matice:= {{q11,q12,q13,q13k},{q12,q22,q23,q23k},{q13,q23,q33,q34},{q13k,q23k,q34,q33k}}

Solve[D[{w1,1-w1-2*w2,w2,w2}.matice.{{w1},{1-w1-2*w2},{w2},{w2}},w1]==0&& D[{w1,1-w1-2*w2,w2,w2}.matice.{{w1},{1-w1-2*w2},{w2},{w2}},w2]==0 , {w1,w2}]

Sw1[mu_,k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_]:= -(((-4 q12[mu,t,sigma1, sigma2]+2 q13[mu,k1,sigma1,sigma2,t,x0,y0]+2 q13[mu,k2,sigma1,sigma2,t,x0,y0]+4 q22[mu,t,sigma1, sigma2]-2 q23[mu,k1,sigma1,sigma2,t,x0,y0]-2 q23[mu,k2,sigma1,sigma2,t,x0,y0]) (-4 q22[mu,t,sigma1, sigma2]+2 q23[mu,k1,sigma1,sigma2,t,x0,y0]+2 q23[mu,k2,sigma1,sigma2,t,x0,y0])-(2 q12[mu,t,sigma1, sigma2]-2 q22[mu,t,sigma1, sigma2]) (4 q22[mu,t,sigma1, sigma2]-6 q23[mu,k1,sigma1,sigma2,t,x0,y0]-6 q23[mu,k2,sigma1,sigma2,t,x0,y0]-2 (-2 q22[mu,t,sigma1, sigma2]+q23[mu,k1,sigma1,sigma2,t,x0,y0]+q23[mu,k2,sigma1,sigma2,t,x0,y0])+2 q33[mu,k1,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]+4 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]))/((((-4 q12[mu,t,sigma1, sigma2]+2 q13[mu,k1,sigma1,sigma2,t,x0,y0]+2 q13[mu,k2,sigma1,sigma2,t,x0,y0]+4 q22[mu,t,sigma1, sigma2]-2 q23[mu,k1,sigma1,sigma2,t,x0,y0]-2 q23[mu,k2,sigma1,sigma2,t,x0,y0]))^2)-(2 q11[mu,t,sigma1, sigma2]-4 q12[mu,t,sigma1, sigma2]+2 q22[mu,t,sigma1, sigma2]) (4 q22[mu,t,sigma1, sigma2]-6 q23[mu,k1,sigma1,sigma2,t,x0,y0]-6 q23[mu,k2,sigma1,sigma2,t,x0,y0]-2 (-2 q22[mu,t,sigma1, sigma2]+q23[mu,k1,sigma1,sigma2,t,x0,y0]+q23[mu,k2,sigma1,sigma2,t,x0,y0])+2 q33[mu,k1,sigma1,sigma2,t,x0,y0]+2 q33[mu,k2,sigma1,sigma2,t,x0,y0]+4 q34[mu,k1,k2,sigma1,sigma2,t,x0,y0])))

Sw1[0.05,2,3,0.3,0.3,10,1,1]

Sw2[mu_,k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_]:=-((-2 ((q12[mu,t,sigma1, sigma2])^2)+q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0]+q12[mu,t,sigma1, sigma2] q13[mu,k2,sigma1,sigma2,t,x0,y0]+2 q11[mu,t,sigma1, sigma2] q22[mu,t,sigma1, sigma2]-q13[mu,k1,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2]-q13[mu,k2,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2]-q11[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0]+q12[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0]-q11[mu,t,sigma1, sigma2] q23[mu,k2,sigma1,sigma2,t,x0,y0]+q12[mu,t,sigma1, sigma2]q23[mu,k2,sigma1,sigma2,t,x0,y0])/(4 ((q12[mu,t,sigma1, sigma2])^2)-4 q12[mu,t,sigma1, sigma2] q13[mu,k1,sigma1,sigma2,t,x0,y0]+((q13[mu,k1,sigma1,sigma2,t,x0,y0])^2)-4 q12[mu,t,sigma1, sigma2] q13[mu,k2,sigma1,sigma2,t,x0,y0]+2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q13[mu,k2,sigma1,sigma2,t,x0,y0]+((q13[mu,k2,sigma1,sigma2,t,x0,y0])^2)-4 q11[mu,t,sigma1, sigma2] q22[mu,t,sigma1, sigma2]+4 q13[mu,k1,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2]+4 q13[mu,k2,sigma1,sigma2,t,x0,y0] q22[mu,t,sigma1, sigma2]+4 q11[mu,t,sigma1, sigma2]q23[mu,k1,sigma1,sigma2,t,x0,y0]-4 q12[mu,t,sigma1, sigma2] q23[mu,k1,sigma1,sigma2,t,x0,y0]-2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0]-2 q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k1,sigma1,sigma2,t,x0,y0]+((q23[mu,k1,sigma1,sigma2,t,x0,y0])^2)+4 q11[mu,t,sigma1, sigma2] q23[mu,k2,sigma1,sigma2,t,x0,y0]-4 q12[mu,t,sigma1, sigma2] q23[mu,k2,sigma1,sigma2,t,x0,y0]-2 q13[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]-2 q13[mu,k2,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]+2 q23[mu,k1,sigma1,sigma2,t,x0,y0] q23[mu,k2,sigma1,sigma2,t,x0,y0]+((q23[mu,k2,sigma1,sigma2,t,x0,y0])^2)-q11[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0]-q22[mu,t,sigma1, sigma2] q33[mu,k1,sigma1,sigma2,t,x0,y0]-q11[mu,t,sigma1, sigma2] q33[mu,k2,sigma1,sigma2,t,x0,y0]+2 q12[mu,t,sigma1, sigma2] q33[mu,k2,sigma1,sigma2,t,x0,y0]-q22[mu,t,sigma1, sigma2] q33[mu,k2,sigma1,sigma2,t,x0,y0]-2 q11[mu,t,sigma1, sigma2] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]+4 q12[mu,t,sigma1, sigma2]q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2 q22[mu,t,sigma1, sigma2] q34[mu,k1,k2,sigma1,sigma2,t,x0,y0]))

Sw2[0.05,2,3,0.3,0.3,10,1,1]

TunedApMaxSW1[mu_,sigma1_,sigma2_,t_,x0_,y0_]:= FindMinimum[{Norm[{Sw1[mu,k1,k2,sigma1,sigma2,t,x0,y0],1-Sw1[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2*Sw2[mu,k1,k2,sigma1,sigma2,t,x0,y0],Sw2[mu,k1,k2,sigma1,sigma2,t,x0,y0],Sw2[mu,k1,k2,sigma1,sigma2,t,x0,y0]}.matka1[mu,k1,k2,sigma1,sigma2,t,x0,y0].{{Sw1[mu,k1,k2,sigma1,sigma2,t,x0,y0]},{1-Sw1[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2*Sw2[mu,k1,k2,sigma1,sigma2,t,x0,y0]},{Sw2[mu,k1,k2,sigma1,sigma2,t,x0,y0]},{Sw2[mu,k1,k2,sigma1,sigma2,t,x0,y0]}}],0<k1<k2},{k1,k2}]

TunedPayoffSW[mu_,sigma1_,sigma2_,k1_,k2_,t_,x0_,y0_]:=NIntegrate[Log[(Sw1[mu,k1,k2,sigma1,sigma2,t,x0,y0]*qx[x,sigma1,t] + (1-Sw1[mu,k1,k2,sigma1,sigma2,t,x0,y0]-2*Sw2[mu,k1,k2,sigma1,sigma2,t,x0,y0])*qy[x,sigma1,t]+Sw2[mu,k1,k2,sigma1,sigma2,t,x0,y0]*pmax1[x,k1,sigma1,t,x0,y0] + Sw2[mu,k1,k2,sigma1,sigma2,t,x0,y0]*pmax1[x,k2,sigma1,t,x0,y0])/qy[x,sigma1,t]]*qm[x,mu,sigma2,t],{x,-Infinity,Infinity}]

TunedPayoffSW[0.05,0.2,0.3,0.76,1.26,1,1,1]

TunedPayoffSW[0.05,0.3,0.3,0.79,1.2727,1,1,1]

TunedApMaxSW1[0.05,0.3,0.3,10,1,1]

TunedApMaxSW1[-0.05,0.3,0.3,10,1,1]

TunedApMaxSW1[0.05,0.2,0.3,1,1,1]

Sw1[0.05,0.76,1.26,0.2,0.3,1,1,1]

Sw2[0.05,0.76,1.26,0.2,0.3,1,1,1]

TunedApMaxSW1[0.05,0.2,0.2,1,1,1]

Densities comparison

Define approximate densities:

TunedDensity0[x_,mu_,sigma1_,sigma2_,t_]:=(qx[x,sigma1,t]*Tw[sigma1,sigma2,t,mu]+(1-Tw[sigma1,sigma2,t,mu])*qy[x,sigma1,t])

TunedDensity1[x_,mu_,k_,sigma1_,sigma2_,t_,x0_,y0_]:= (Tw1[mu,k,sigma1,sigma2,t,x0,y0]*qx[x,sigma1,t] + Tw2[mu,k,sigma1,sigma2,t,x0,y0]qy[x,sigma1,t] + (1-Tw1[mu,k,sigma1,sigma2,t,x0,y0]-Tw2[mu,k,sigma1,sigma2,t,x0,y0])*pmax1[x,k,sigma1,t,x0,y0])

TunedDensity2[x_,mu_,k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_]:= (v1[mu,k1,k2,sigma1,sigma2,t,x0,y0]*qx[x,sigma1,t] + v2[mu,k1,k2,sigma1,sigma2,t,x0,y0]qy[x,sigma1,t] + v3[mu,k1,k2,sigma1,sigma2,t,x0,y0]*pmax1[x,k1,sigma1,t,x0,y0]+ (1-v1[mu,k1,k2,sigma1,sigma2,t,x0,y0]-v2[mu,k1,k2,sigma1,sigma2,t,x0,y0]-v3[mu,k1,k2,sigma1,sigma2,t,x0,y0])*pmax1[x,k2,sigma1,t,x0,y0])

MarkowitzDensity0[x_,mu_,sigma1_,sigma2_,t_]:=vaha[mu,sigma1,sigma2,t]*px[x,sigma1,t] + (1-vaha[mu,sigma1,sigma2,t])*py[x,sigma1,t]

MarkowitzDensity1[x_,k_,mu_,sigma1_,sigma2_,t_,x0_,y0_]:= (w1[k,sigma1,sigma2,t,x0,y0,mu]*px[x,sigma1,t] + w2[k,sigma1,sigma2,t,x0,y0,mu]py[x,sigma1,t] + (1-w1[k,sigma1,sigma2,t,x0,y0,mu]-w2[k,sigma1,sigma2,t,x0,y0,mu])*pmax[x,k,sigma1,t,x0,y0])

MarkowitzDensity2[x_,k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_,mu_]:= (Mw1[k1,k2,sigma1,sigma2,t,x0,y0,mu]*px[x,sigma1,t] + Mw2[k1,k2,sigma1,sigma2,t,x0,y0,mu]*py[x,sigma1,t] + Mw3[k1,k2,sigma1,sigma2,t,x0,y0,mu]*pmax[x,k1,sigma1,t,x0,y0] +Mw4[k1,k2,sigma1,sigma2,t,x0,y0,mu]*pmax[x,k2,sigma1,t,x0,y0] )

HardDensity0[x_,mu_,sigma1_,sigma2_,t_,w_]:=(w*qx[x,sigma1,t ]+ (1-w)*qy[x,sigma1,t])

HardDensity1[x_,mu_,k_,sigma1_,sigma2_,t_,x0_,y0_,w1_,w2_]:=(w1*qx[x,sigma1,t] + w2*qy[x,sigma1,t] + (1-w1-w2)*pmax1[x,k,sigma1,t,x0,y0])

HardDensity2[x_,mu_,k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_,w1_,w2_,w3_]:=(w1*qx[x,sigma1,t ]+ w2*qy[x,sigma1,t] +  w3*pmax1[x,k1,sigma1,t,x0,y0]+ (1-w1-w2-w3)*pmax1[x,k2,sigma1,t,x0,y0])

Plotting them for all three cases:

Density0Fit[mu_,sigma1_,sigma2_,t_,w_]:= Plot[{MarkowitzDensity0[x,mu,sigma1,sigma2,t],HardDensity0[x,mu,sigma1,sigma2,t,w],TunedDensity0[x,mu,sigma1,sigma2,t], qm[x,mu,sigma2,t]}, {x,-5,5},PlotLabel ->"No contract portfolios and market density comparison",PlotLegends->{"Markowitz density","Numerical density","Tuned density", "market opinion density \!\(\*TemplateBox[<|"boxes" -> FormBox[
SubscriptBox[
StyleBox["p", "TI"], 
StyleBox["M", "TI"]], TraditionalForm], "errors" -> {}, "input" -> "p_M", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)"},PlotRange->Automatic]

Density0Fit[0.02,0.2,0.2,10,0.767]

Density1Fit[mu_,kH_,kT_,kM_,sigma1_,sigma2_,t_,x0_,y0_,w1_,w2_]:=Plot[{MarkowitzDensity1[x,kM,mu,sigma1,sigma2,t,x0,y0],HardDensity1[x,mu,kH,sigma1,sigma2,t,x0,y0,w1,w2],TunedDensity1[x,mu,kT,sigma1,sigma2,t,x0,y0], qm[x,mu,sigma2,t]}, {x,-5,5},PlotLabel ->"One contract ortfolios and market density comparison",PlotLegends->{"Markowitz density","Numerical density","Tuned density", "market opinion density \!\(\*TemplateBox[<|"boxes" -> FormBox[
SubscriptBox[
StyleBox["p", "TI"], 
StyleBox["M", "TI"]], TraditionalForm], "errors" -> {}, "input" -> "p_M", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)"},PlotRange->Automatic]

Density1Fit[0.05,1.023,1.019,1.681,0.3,0.3,10,1,1,0.937,0.84]

Density2Fit[mu_,kH1_,kH2_,kT1_,kT2_,kM1_,kM2_,sigma1_,sigma2_,t_,x0_,y0_,w1_,w2_,w3_]:= Plot[{MarkowitzDensity2[x,kM1,kM2,sigma1,sigma2,t,x0,y0,mu],HardDensity2[x,mu,kH1,kH2,sigma1,sigma2,t,x0,y0,w1,w2,w3],TunedDensity2[x,mu,kT1,kT2,sigma1,sigma2,t,x0,y0], qm[x,mu,sigma2,t]}, {x,-5,5},PlotLabel ->"Two contracts portfolios and market density comparison",PlotLegends->{"Markowitz density","Numerical density","Tuned density", "market opinion density \!\(\*TemplateBox[<|"boxes" -> FormBox[
SubscriptBox[
StyleBox["p", "TI"], 
StyleBox["M", "TI"]], TraditionalForm], "errors" -> {}, "input" -> "p_M", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)"},PlotRange->Automatic]

Density2Fit[0.05,2.411,4.498,0.279,4.087,0.1,1.328,0.3,0.4,10,1,1,0.294,-1.477,2.184]

Relative entropy in all cases

NO CONTRACT

DKLH0[mu_,sigma1_,sigma2_,t_,w_]:=NIntegrate[qm[x,mu,sigma2,t]*Log[qm[x,mu,sigma2,t]/HardDensity0[x,mu,sigma1,sigma2,t,w]],{x,-Infinity,Infinity}]

DKLH0[0.03,0.2,0.2,10,0.767]

DKLM0[mu_,sigma1_,sigma2_,t_]:=NIntegrate[qm[x,mu,sigma2,t]*Log[qm[x,mu,sigma2,t]/MarkowitzDensity0[x,mu,sigma1,sigma2,t]],{x,-Infinity,Infinity}]

DKLM0[0.03,0.2,0.2,10]

DKLT0[mu_,sigma1_,sigma2_,t_]:=NIntegrate[qm[x,mu,sigma2,t]*Log[qm[x,mu,sigma2,t]/TunedDensity0[x,mu,sigma1,sigma2,t]],{x,-Infinity,Infinity}]

DKLT0[0.03,0.2,0.2,10]

ONE CONTRACT

DKLH1[mu_,sigma1_,sigma2_,t_,k_,x0_,y0_,w1_,w2_]:= NIntegrate[qm[x,mu,sigma2,t]*Log[qm[x,mu,sigma2,t]/HardDensity1[x,mu,k,sigma1,sigma2,t,x0,y0,w1,w2]],{x,-Infinity,Infinity}]

DKLH1[0.05,0.3,0.3,10,1.023,1,1,0.937,0.84]

DKLM1[mu_,sigma1_,sigma2_,t_,k_,x0_,y0_]:= NIntegrate[qm[x,mu,sigma2,t]*Log[qm[x,mu,sigma2,t]/MarkowitzDensity1[x,k,mu,sigma1,sigma2,t,x0,y0]],{x,-Infinity,Infinity}]

DKLM1[0.05,0.3,0.3,10,1.681,1,1]

DKLT1[mu_,sigma1_,sigma2_,t_,k_,x0_,y0_]:= NIntegrate[qm[x,mu,sigma2,t]*Log[qm[x,mu,sigma2,t]/TunedDensity1[x,mu,sigma1,sigma2,t,k,x0,y0]],{x,-Infinity,Infinity}]

DKLT1[0.05,0.3,0.3,10,1.023,1,1]

TWO CONTRACTS

DKLH2[mu_,k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_,w1_,w2_,w3_] := NIntegrate[qm[x,mu,sigma2,t]*Log[qm[x,mu,sigma2,t]/HardDensity2[x,mu,k1,k2,sigma1,sigma2,t,x0,y0,w1,w2,w3]],{x,-Infinity,Infinity}]

DKLH2[0.05,2.411,2.875,0.3,0.4,10,1,1,0.293,-1.477,-0.059]

DKLM2[mu_,k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_] := NIntegrate[qm[x,mu,sigma2,t]*Log[qm[x,mu,sigma2,t]/MarkowitzDensity2[x,k1,k2,sigma1,sigma2,t,x0,y0,mu]],{x,-Infinity,Infinity}]

DKLM2[0.05,0.1,1.328,0.3,0.4,10,1,1]

DKLT2[mu_,k1_,k2_,sigma1_,sigma2_,t_,x0_,y0_] := NIntegrate[qm[x,mu,sigma2,t]*Log[qm[x,mu,sigma2,t]/TunedDensity2[x,mu,k1,k2,sigma1,sigma2,t,x0,y0]],{x,-Infinity,Infinity}]

DKLT2[0.05,0.279,4.087,0.3,0.4,10,1,1]

completeDKL[mu_,sigma1_,sigma2_,t_] := NIntegrate[qm[x,mu,sigma2,t]*Log[qm[x,mu,sigma2,t]/qy[x,sigma1,t]],{x,-Infinity,Infinity}]

completeDKL[0.03,0.2,0.3,50]

Option strategies

Define prices an P/L functions and plot them:

CallOptionPrice[sigma1_,t_,xy_,k_]:= xy*CDF[NormalDistribution[0,1], sigma1*Sqrt[t]/2+1/(sigma1*Sqrt[t])*Log[xy/k]] - k*CDF[NormalDistribution[0,1], -sigma1*Sqrt[t]/2+1/(sigma1*Sqrt[t])*Log[xy/k]]

PutOptionPrice[sigma1_,t_,xy_,k_]:= CallOptionPrice[sigma1,t,xy,k]+k-xy

BullishSpreadPayoff[p_,sigma1_,t_,xy_,k1_,k2_]:= -CallOptionPrice[sigma1,t,xy,k1]*Boole[p<k1] + Boole[p>=k1]*(p-CallOptionPrice[sigma1,t,xy,k1]-k1) -( Boole[p>=k2]*(p-k2-CallOptionPrice[sigma1,t,xy,k2]) -CallOptionPrice[sigma1,t,xy,k2]* Boole[p<k2])

option[p_,sigma1_,t_,xy_,k2_]:=-( Boole[p>=k2]*(p-k2-CallOptionPrice[sigma1,t,xy,k2]) -CallOptionPrice[sigma1,t,xy,k2]* Boole[p<k2])

Plot[BullishSpreadPayoff[p,0.3,1,1,0.46,2.11],{p,0,2.5},PlotLabel->"Payoff of Bullish Spread with optimal strikes",AxesLabel->{"price", "P/L"}]

Int[sigma1_,t_,xy_,k1_,k2_] :=NIntegrate[Max[xy*qx[x,sigma1,t] - k1*qy[x,sigma1,t],0] + Min[-xy*qx[x,sigma1,t] + k2*qy[x,sigma1,t],0] + (CallOptionPrice[sigma1,t,xy,k1]-CallOptionPrice[sigma1,t,xy,k2])*qy[x,sigma1,t],{x,-Infinity,Infinity}]

DensityBullishSpread[x_,sigma1_,t_,xy_,k1_,k2_] := (Max[xy*qx[x,sigma1,t] - k1*qy[x,sigma1,t],0] + Min[-xy*qx[x,sigma1,t] + k2*qy[x,sigma1,t],0] + (CallOptionPrice[sigma1,t,xy,k1]-CallOptionPrice[sigma1,t,xy,k2])*qy[x,sigma1,t])

Plot[{DensityBullishSpread[x,0.3,10,1,0.46,2.11]/0.90269,qx[x,0.3,10],qy[x,0.3,10]},{x,-3,3.5},PlotLegends->{"Bullish spread density", "SPD of X","SPD of Y"}, PlotLabel->"SPD of Bullish spread"]

BearishSpreadpayoff[p_,sigma1_,t_,xy_,k2_,k1_]:= PutOptionPrice[sigma1,t,xy,k2]*Boole[p>=k2] + Boole[p<k2]*(p+PutOptionPrice[sigma1,t,xy,k2]-k2) +( Boole[p<k1]*(-p+k1-PutOptionPrice[sigma1,t,xy,k1]) -PutOptionPrice[sigma1,t,xy,k1]* Boole[p>=k1])

Putoption[p_,sigma1_,t_,xy_,k1_]:= ( Boole[p<k1]*(-p+k1-PutOptionPrice[sigma1,t,xy,k1]) -PutOptionPrice[sigma1,t,xy,k1]* Boole[p>=k1])

Plot[Putoption[p,0.3,10,1,0.677],{p,0,2.5},PlotRange->{-0.3,0.6},PlotLabel->"P/L of a put option with optimal strike",AxesLabel->{"price","P/L"}]

Plot[BearishSpreadpayoff[p,0.3,10,1,0.34,0.84],{p,0,1.8},PlotLabel->"Payoff of Bearish Spread with optimal strikes",AxesLabel->{"price","P/L"}]

LongStranglepayoff[p_,sigma1_,t_,xy_,k1_,k2_]:= -PutOptionPrice[sigma1,t,xy,k1]*Boole[p>=k1] + Boole[p<k1]*(-p-PutOptionPrice[sigma1,t,xy,k1]+k1) +( Boole[p>=k2]*(p-k2-CallOptionPrice[sigma1,t,xy,k2]) -CallOptionPrice[sigma1,t,xy,k2]* Boole[p<k2])

option[p_,sigma1_,t_,xy_,k1_]:=-PutOptionPrice[sigma1,t,xy,k1]*Boole[p>=k1] + Boole[p<k1]*(-p-PutOptionPrice[sigma1,t,xy,k1]+k1)

option[p_,sigma1_,t_,xy_,k1_]:= Boole[p<k1]*(p-k1-CallOptionPrice[sigma1,t,xy,k1]) -CallOptionPrice[sigma1,t,xy,k1]* Boole[p<k1]

Plot[LongStranglepayoff[p,0.3,1,1,1.0014,1.56493],{p,0,2.5},PlotLabel->"Payoff of Long Strangle with optimal strikes"]

Weights comparison

We print out all the weight in one place for convenience.

No contract

same volatility

HardMaximum[0.03,0.2,0.2,10]

vaha[0.03,0.2,0.2,10]

MarkowitzPayoff0[0.03,0.2,0.2,10]

Tw[0.2,0.2,10,0.03]

TunedPayoff0[0.2,0.2,10,0.03]

Here we plot out the payoff comparison in time from Section 3.1

Plot[{TunedPayoff0[0.2,0.2,t,0.03],MarkowitzPayoff0[0.03,0.2,0.2,t],(0.03/0.2)^2*t/2},{t,0,100},PlotLegends->{"Tuned payoff", "Markowitz payoff","Optimal payoff"},PlotLabel->"Payoff comparison"]

one contract

HardMaximum1[0.05,0.3,0.3,10,1,1]

MarkowitzApMax1[0.05,0.3,0.3,10,1,1]

w1[1.68055,0.3,0.3,10,1,1,0.05]

w2[1.68055,0.3,0.3,10,1,1,0.05]

MarkowitzPayoff1[1.68055,0.05,0.3,0.3,10,1,1]

TunedApMax1[0.05,0.3,0.3,10,1,1]

Tw1[0.05,1.019,0.3,0.3,10,1,1]

Tw2[0.05,1.019,0.3,0.3,10,1,1]

TunedPayoff1[0.05,0.3,0.3,10,1.019,1,1]

two contracts portfolio

HardMaximum2[0.05,0.3,0.4,10,1,1]

HardPayoff2[mu_,sigma1_,sigma2_,t_,x0_,y0_, k1_,k2_,w1_,w2_,w3_]:= NIntegrate[qm[x,mu,sigma2,t]*Log[(w1*qx[x,sigma1,t ]+ w2*qy[x,sigma1,t] +  w3*pmax1[x,k1,sigma1,t,x0,y0]+ (1-w1-w2-w3)*pmax1[x,k2,sigma1,t,x0,y0])/qy[x,sigma1,t]],{x,-Infinity,Infinity}]

TunedApMax2[0.05,0.3,0.4,10,1,1]

v1[0.05,0.279,4.087,0.3,0.4,10,1,1]

v2[0.05,0.279,4.087,0.3,0.4,10,1,1]

v3[0.05,0.279,4.087,0.3,0.4,10,1,1]

TunedPayoff2[0.05,0.279,4.087,0.3,0.4,10,1,1]

MarkowitzApMax2[0.3,0.4,10,1,1,0.05]

Mw1[0.1,1.388,0.3,0.4,10,1,1,0.05]

Mw2[0.1,1.388,0.3,0.4,10,1,1,0.05]

Mw3[0.1,1.388,0.3,0.4,10,1,1,0.05]

MarkowitzPayoff2[0.1,1.388,0.3,0.4,10,1,1,0.05]

Mean variance plots

points

pointH[mu_,k_,sigma1_,sigma2_,t_,x0_,y0_,w1_,w2_] := {Sqrt[NIntegrate[(HardDensity1[x,mu,k,sigma1,sigma2,t,x0,y0,w1,w2]/qy[x,sigma1,t]-1)^2*qy[x,sigma1,t],{x,-Infinity,Infinity}]],NIntegrate[(HardDensity1[x,mu,k,sigma1,sigma2,t,x0,y0,w1,w2]/qy[x,sigma1,t]-1)*qm[x,mu,sigma2,t],{x,-Infinity,Infinity}] }

pointH[0.05,1.023,0.3,0.3,10,1,1,0.937,0.84]

pointH[0.05,1.02,0.3,0.35,1,1,1,-0.5133,-0.75224]

pointT[mu_,k_,sigma1_,sigma2_,t_,x0_,y0_] := {Sqrt[NIntegrate[(TunedDensity1[x,mu,k,sigma1,sigma2,t,x0,y0]/qy[x,sigma1,t]-1)^2*qy[x,sigma1,t],{x,-Infinity,Infinity}]],NIntegrate[(TunedDensity1[x,mu,k,sigma1,sigma2,t,x0,y0]/qy[x,sigma1,t]-1)*qm[x,mu,sigma2,t],{x,-Infinity,Infinity}] }

pointT[0.05,1.019,0.3,0.3,10,1,1]

pointT[0.05,1.017,0.3,0.35,1,1,1]

pointMark[mu_,k_,sigma1_,sigma2_,t_,x0_,y0_] := {Sqrt[NIntegrate[(MarkowitzDensity1[x,k,mu,sigma1,sigma2,t,x0,y0]/qy[x,sigma1,t]-1)^2*qy[x,sigma1,t],{x,-Infinity,Infinity}]],NIntegrate[(MarkowitzDensity1[x,k,mu,sigma1,sigma2,t,x0,y0]/qy[x,sigma1,t]-1)*qm[x,mu,sigma2,t],{x,-Infinity,Infinity}] }

pointMark[0.05,1.681,0.3,0.3,10,1,1]

pointMark[0.05,1.2485,0.3,0.35,1,1,1]

pointM[mu_,sigma1_,sigma2_,t_,x0_,y0_] := {Sqrt[NIntegrate[(qm[x,mu,sigma2,t]/qy[x,sigma1,t]-1)^2*qy[x,sigma1,t],{x,-Infinity,Infinity}]],NIntegrate[(qm[x,mu,sigma2,t]/qy[x,sigma1,t]-1)*qm[x,mu,sigma2,t],{x,-Infinity,Infinity}] }

pointM[0.05,0.3,0.3,10,1,1]

pointM[0.05,0.3,0.35,1,1,1]

pointX[mu_,sigma1_,sigma2_,t_,x0_,y0_] := {Sqrt[NIntegrate[(qx[x,sigma1,t]/qy[x,sigma1,t]-1)^2*qy[x,sigma1,t],{x,-Infinity,Infinity}]],NIntegrate[(qx[x,sigma1,t]/qy[x,sigma1,t]-1)*qm[x,mu,sigma2,t],{x,-Infinity,Infinity}] }

pointX[0.05,0.3,0.3,10,1,1]

pointX[0.05,0.3,0.35,1,1,1]

frontiers:

maxFrontier[k_,mu_,sigma1_,sigma2_,t_,x0_,y0_]:= {Sqrt[NIntegrate[(pmax1[x,k,sigma1,t,x0,y0]/qy[x,sigma1,t]-1)^2*qy[x,sigma1,t],{x,-Infinity,Infinity}]],NIntegrate[(pmax1[x,k,sigma1,t,x0,y0]/qy[x,sigma1,t]-1)*qm[x,mu,sigma2,t],{x,-Infinity,Infinity}] };

maxFrontierK[k_]:=maxFrontier[k,0.05,0.3,0.3,10,1,1]

maxFrontierK2[k_]:=maxFrontier[k,0.05,0.3,0.35,1,1,1]

maxFrontier[1,0.05,0.3,0.3,10,1,1]

VDMarkowitz1[k_,mu_,sigma1_,sigma2_,t_,x0_,y0_]:={Sqrt[NIntegrate[qy[x,sigma1,t]*((MarkowitzDensity1[x,k,mu,sigma1,sigma2,t,x0,y0])/qy[x,sigma1,t]-1)^2,{x,-Infinity,Infinity}]],NIntegrate[qm[x,mu,sigma2,t]*((MarkowitzDensity1[x,k,mu,sigma1,sigma2,t,x0,y0])/qy[x,sigma1,t]-1),{x,-Infinity,Infinity}]}

Markowitz1[k_]:= VDMarkowitz1[k,0.05,0.3,0.3,10,1,1]

Markowitz2[k_]:= VDMarkowitz1[k,0.05,0.3,0.35,1,1,1]

xmaxfrontier[w_] := {Sqrt[Dot[{w,0,1-w},B[1.023,0.3,0.3,10,1,1,0.05].{w,0,1-w}]],w*mu1[0.3,0.3,0.05,10] + (1-w)*mu3[1.023,0.3,0.3,10,1,1,0.05]}

xmaxfrontier2[w_] := {Sqrt[Dot[{w,0,1-w},B[1.02,0.3,0.35,1,1,1,0.05].{w,0,1-w}]],w*mu1[0.3,0.35,0.05,1] + (1-w)*mu3[1.02,0.3,0.35,1,1,1,0.05]}

ListLinePlot[Array[xmaxfrontier,30,{-4,6}],AxesLabel->{"Volatility","Drift"},PlotRange->Automatic]

B[1.681,0.3,0.3,10,1,1,0.05]

{{1.4596,0.520761},{0.520761,0.221898}}

B[1.023,0.3,0.3,10,1,1,0.05]

{{1.4596`,0.879233},{0.879233,0.56549}}

B[1.2485,0.3,0.35,1,1,1,0.05]

{{0.0941742,0.024141},{0.024141,0.01088}}

parabaolaM[x_]:= {x,x*x/2 + Dot[{mu1[0.3,0.3,0.05,10],mu3[1.681,0.3,0.3,10,1,1,0.05]},Inverse[{{1.459603,0.5207613},{0.5207613,0.221898}}].{mu1[0.3,0.3,0.05,10],mu3[1.681,0.3,0.3,10,1,1,0.05]}]}

parabaolaM2[x_]:= {x,x*x/2 + 1/2*Dot[{mu1[0.3,0.35,0.05,1],mu3[1.2485,0.3,0.35,1,1,1,0.05]},Inverse[{{0.0941742,0.0241413},{0.0241413`,0.01088}}].{mu1[0.3,0.35,0.05,1],mu3[1.2485,0.3,0.35,1,1,1,0.05]}]}

parabaolaH[x_]:= {x,x*x/2 + Dot[{mu1[0.3,0.3,0.05,10],mu3[1.023,0.3,0.3,10,1,1,0.05]},Inverse[{{1.45960,0.8792330},{0.8792330,0.56549}}].{mu1[0.3,0.3,0.05,10],mu3[1.023,0.3,0.3,10,1,1,0.05]}]}

parabaolaH2[x_]:= {x,x*x/2 +(0.084839288151-1/2*0.27581*0.27581) }

We plot the case from Section 3.2:

ListLinePlot[{Array[maxFrontierK, 100,{0.05,13}],
{{0,0},4*pointH[0.05,1.023,0.3,0.3,10,1,1,0.937,0.84]},Array[xmaxfrontier,100,{-4,5}],Array[parabaolaM,100,{0,1}]},
Frame->True,FrameLabel->{"Portfolio Volatility","Portfolio Return"},PlotLabel->"Mean-Variance Frontiers",PlotLegends->{"Contracts","Tangent","Without Y","Markowitz app. logreturn"}
,Epilog->{Point[pointH[0.05,1.023,0.3,0.3,10,1,1,0.937,0.84]],Text["H",{0.59,0.36}]
,Point[pointT[0.05,1.019,0.3,0.3,10,1,1]],Text["M",{0.57,0.33}]
,Point[pointMark[0.05,1.681,0.3,0.3,10,1,1]],Text["V",{0.558,0.295}],Point[pointM[0.05,0.3,0.3,10,1,1]],Text["F",{0.58,0.31}],Point[maxFrontierK[1.023]],Text["MH",{0.75,0.38}] },PlotRange->{{0.34,0.85},{0,0.5}}]

We plot the more transparent case:

ListLinePlot[{Array[maxFrontierK2, 100,{0.05,13}],
{{0,0},2*pointH[0.05,1.02,0.3,0.35,1,1,1,-0.5133,-0.75224]},Array[xmaxfrontier2,100,{-4,5}],Array[parabaolaM2,100,{0,1}],Array[parabaolaH2,100,{0,1}]},
Frame->True,FrameLabel->{"Portfolio Volatility","Portfolio Return"},PlotLabel->"Mean-Variance Frontiers",PlotLegends->{(*"2 Assets",*)"Contracts","Tangent","Without Y","Markowitz app. logreturn"},Epilog->{(*Points labeled as 1,2,and 3*)Point[pointH[0.05,1.02,0.3,0.35,1,1,1,-0.5133,-0.75224]],Text["H",{0.28,0.1}]
,Point[pointT[0.05,1.017,0.3,0.35,1,1,1]],Text["F",{0.27,0.09}]
,Point[pointMark[0.05,1.2485,0.3,0.35,1,1,1]],Text["V",{0.32,0.105}],Point[pointM[0.05,0.3,0.35,1,1,1]],Text["M",{0.355,0.125}],Point[maxFrontierK2[1.02]],Text["MH",{0.185,0.045}],Point[pointX[0.05,0.3,0.35,1,1,1]],Text["X",{0.31,0.075}],Point[{0,0}],Text["Y",{0.01,0.01}] },PlotRange->{{0,0.4},{0,0.15}}
]

2 Asset payoff optimum comparison

We plot the payoffs of all three methods in a 2 asset portfolio as functions of weight:

vynosH[w_,mu_,sigma1_,sigma2_,t_]:=NIntegrate[qm[x,mu,sigma2,t]*Log[(w*qx[x,sigma1,t ]+ (1-w)*qy[x,sigma1,t])/qy[x,sigma1,t]],{x,-Infinity,Infinity}]

vynosM[w_,mu_,sigma1_,sigma2_,t_]:= w*mu1[sigma1,sigma2,mu,t]-1/2*{w,1-w}.{{-1+(E^((sigma1^2) t)),0},{0,0}}.{{w},{1-w}}

vynosT[w_,mu_,sigma1_,sigma2_,t_]:= 1/2*(mu/sigma2)^2*t-1/2{w,1-w}.{{q11[mu,t,sigma1, sigma2],q12[mu,t,sigma1, sigma2]},{q12[mu,t,sigma1, sigma2],q22[mu,t,sigma1, sigma2]}}.{{w},{1-w}}

graf = Plot[{vynosH[w,0.03,0.2,0.2,10],vynosM[w,0.03,0.2,0.2,10],vynosT[w,0.03,0.2,0.2,10]},{w,0,1.1},PlotLegends ->{"Exact payoff","Markowitz payoff","Fischer Payoff"},PlotLabel->"The axact and approximate payoffs with their solutions",AxesLabel->{"weight",""}]

We add the optimal potfolios obtained by every method to each curve:

body = {{0.767,0.11},{vaha[0.03,0.2,0.2,10],Norm[vynosM[vaha[0.03,0.2,0.2,10],0.03,0.2,0.2,10]]},{Tw[0.2,0.2,10,0.03],Norm[vynosT[Tw[0.2,0.2,10,0.03],0.03,0.2,0.2,10]]}}

body = {{0.767,0.11},{vaha[0.03,0.2,0.2,10],Norm[vynosM[vaha[0.03,0.2,0.2,10],0.03,0.2,0.2,10]]},{Tw[0.2,0.2,10,0.03],Norm[vynosT[Tw[0.2,0.2,10,0.03],0.03,0.2,0.2,10]]},{0.767,Norm[vynosM[0.767,0.03,0.2,0.2,10]]},{Tw[0.2,0.2,10,0.03],Norm[vynosM[Tw[0.2,0.2,10,0.03],0.03,0.2,0.2,10]]}, {vaha[0.03,0.2,0.2,10],Norm[vynosH[vaha[0.03,0.2,0.2,10],0.03,0.2,0.2,10]]},{Tw[0.2,0.2,10,0.03],Norm[vynosH[Tw[0.2,0.2,10,0.03],0.03,0.2,0.2,10]]},{0.767,Norm[vynosT[0.767,0.03,0.2,0.2,10]]},{vaha[0.03,0.2,0.2,10],Norm[vynosT[vaha[0.03,0.2,0.2,10],0.03,0.2,0.2,10]]}};
colors = {Blue,Orange,Green,Blue,Green,Orange,Green,Blue,Orange};
labeledData=MapThread[Style[#1,#2]&,{body,colors}];
graf1 = ListPlot[labeledData,PlotStyle->PointSize[0.01]]

Show[graf,graf1]

graf = Plot[{vynosH[w,0.05,0.3,0.3,10],vynosM[w,0.05,0.3,0.3,10],vynosT[w,0.05,0.3,0.3,10]},{w,0,1.1},PlotLegends ->{"Exact payoff","Markowitz payoff","Fischer Payoff"},PlotLabel->"The axact and approximate payoffs"]

body = {{0.565,0.129},{vaha[0.05,0.3,0.3,10],Norm[vynosM[vaha[0.05,0.3,0.3,10],0.05,0.3,0.3,10]]},{Tw[0.3,0.3,10,0.05],Norm[vynosT[Tw[0.3,0.3,10,0.05],0.05,0.3,0.3,10]]}};
nazvy = {"H, F","M",""};
labeledData=MapThread[Labeled[#1,#2]&,{body,nazvy}];
graf1 = ListPlot[labeledData]

Show[graf,graf1]