(* Code to compute optimal indifference prices and optimal position sizes for a put option in the basis risk model*)
(* dS_t / S_t = μ dt + σ(ρdW_t+ (((1-(ρ^2)))^(1/2))dB_t);  dY_t / Y_t =b dt + adW_t
h(y) = (SuperPlus[(Stk-Y_T)])*)
(*U(x) = -(1/α)Exp[-αx]*)
(*M = numerical integration bound for the normal pdf*)
(*Tol = error tolerance for finding the fixed point*)

Stk=20;
T=1;
b=.05;
μ=.06;
a=.3;
σ=.4;
y=15;
M=20;
Tol=.00001;
α=3;

(*GetBetaVal takes in the correlation ρ, the purchase price Px, and tries to solve the equation*)
	(* Px = ((E[h(Y_T)Exp[-γh(Y_T)]])/(E[Exp[-γh(Y_T)]])) where dY_t/Y_t = (b-ρaμ/σ)dt + adW_t*)
(*for γ. Since the right hand side above is monotonic in γ, the optimal (SuperStar[γ]) is found using a simple bifurcation method taking StartHighGamma and StartLowGamma as initial inputs.*)
(*Once the optimal (SuperStar[γ]) is found, the optimal purchase quantitiy (SuperStar[q])=(SuperStar[γ])/(α(1-(ρ^2))) as well as the indifference price p((SuperStar[q])) = -(1/(SuperStar[γ]))Log E[Exp[-(SuperStar[γ]) h(Y_T)]] are computed.*)
(*The outputs are:
   OV1 : optimal (SuperStar[γ]);
	 OV2 : indifference price;
	 OV3 : optimal purchase quantity;
	 OV4 : first order approx. to the indifference price;
     Px : purchase price - also the marginal price for the optimal quantitity;
     expec: E[h(Y_T)] (unique arb free price when ρ = 1);
     counter : counts the iterations until convergence.
    *)

GetGammaVal[ρ_,StartHighGamma_,StartLowGamma_,Px_]:=Module[{hatz,Counter,Error,TempHighGamma,TempLowGamma,TempGamma,TempPx,expec, var,OV1,OV2,OV3,OV4},
Clear[hatz,Counter,Error,TempHighGamma,TempLowGamma,TempGamma,TempPx,expec,var,OV1,OV2,OV3,OV4];
hatz=(1/(a×(T^(1/2))))×(Log[(Stk/y)]-(b-ρ×μ×(a/σ)-((a^2)/2))×T);
f[γ_]:=NIntegrate[(Stk-y×Exp[(b-ρ×μ×(a/σ)-((a^2)/2))×T+a×(T^(1/2))×z])×Exp[-γ×(Stk-y×Exp[(b-ρ×μ×(a/σ)-((a^2)/2))×T+a×(T^(1/2))×z])]×(1/(((2×π))^(1/2)))×Exp[((-(z^2))/2)],{z,-M,hatz}];
g[γ_]:=NIntegrate[Exp[-γ×(Stk-y×Exp[(b-ρ×μ×(a/σ)-((a^2)/2))×T+a×(T^(1/2))×z])]×(1/(((2×π))^(1/2)))×Exp[((-(z^2))/2)],{z,-M,hatz}]+ 1- CDF[NormalDistribution[0,1],hatz];
expec=NIntegrate[(Stk-y×Exp[(b-ρ×μ×(a/σ)-((a^2)/2))×T+a×(T^(1/2))×z])×(1/(((2×π))^(1/2)))×Exp[((-(z^2))/2)],{z,-M,hatz}];
var=NIntegrate[(((Stk-y×Exp[(b-ρ×μ×(a/σ)-((a^2)/2))×T+a×(T^(1/2))×z]))^2)×(1/(((2×π))^(1/2)))×Exp[((-(z^2))/2)],{z,-M,hatz}]-(expec^2);
Counter=0;
Error = 2×Tol;
TempHighGamma=StartHighGamma;
TempLowGamma=StartLowGamma;
While[Error>Tol&&Counter<500,
TempGamma=(TempHighGamma+TempLowGamma)/2;
TempPx=((f[TempGamma])/(g[TempGamma]));
Which[TempPx > Px,TempLowGamma=TempGamma,True,TempHighGamma=TempGamma];
Error = Abs[TempPx-Px];
Counter=Counter+1;
];
Which[Counter==500,OV1="No Convergence";OV2="";OV3="";OV4="",True,OV1=TempGamma//N;OV2=((-1)/TempGamma)×Log[g[TempGamma]];OV3=(TempGamma/(α×(1-(ρ^2))));OV4=expec-((TempGamma/2))×var];
{OV1,OV2,OV3,OV4,Px, expec,Counter}];

RhoMap[ρ_]:=GetGammaVal[ρ,5,-5,5];

DataTable=Table[{.8+j/2000,RhoMap[.8+j/2000]},{j,0,399}];

Needs["PlotLegends`"]

PxTable=ListLinePlot[{Table[{DataTable[[i]][[1]],DataTable[[i]][[2]][[2]]},{i,1,400}],Table[{DataTable[[i]][[1]],DataTable[[400]][[2]][[6]]},{i,1,400}](*, Table[{DataTable[[i]][[1]],DataTable[[i]][[2]][[5]]},{i,1,400}]*)},PlotRange->All,PlotStyle->{Dashing[None],Dashing[Tiny], Dashing[Large]},AxesOrigin->{.8,5},AxesLabel->{"Correlation (Subscript[ρ, n])", "Prices"}(*,PlotLegend->{"(p^n)((OverHat[q]))", "(OverHat[p])", "p"}, LegendShadow->None, LegendBorder-> None,LegendPosition->{1,-.5},LegendSize->{.5,.5}*)]

PosTable=ListLinePlot[{Table[{DataTable[[i]][[1]],DataTable[[i]][[2]][[3]]},{i,1,400}]},PlotStyle->{Dashing[None]},PlotRange->{0,3},AxesOrigin->{.8,0}, AxesLabel->{"Correlation (Subscript[ρ, n])","Position (Subscript[(SuperStar[q]), n])"}]