The T-date equilibrium in an incomplete market

Basak-Cuoco example
Solved on a binomial tree

For non binomial trees, another program should be used.

August 4, 2008
Purely editorail: replaced all instances of "transform" or "Transform" by "Undiscount"

August 4, 2008
Removed function ϕ2F2 from consideration and only used F2

August 4, 2008
Switched to the interpolation of the undiscounted value function

The "jewel" program!
July 25, 2008
Extension to any number of unknowns in preparation for Heaton and Lucas

July 24, 2008
Extension to any number of iterated functions
From July_23_2008_Boston_BD.nb

July 23, 2008
Defined a function ϕ2F2 from the array ϕ2F2T

July 22, 2008
Added some documentation to July_ 14 _ 2008 _Boston.nb

Off[General::spell1]; Off[General::spell];

Needs["PlotLegends`"]

Define state prices (one should not have to open this)

ϕ2[τ_, c2_] := ρ^(T - τ)*c2^(γ2 - 1); 
ϕ1[τ_, c1_] := ρ^(T - τ)*c1^(γ1 - 1);

t = T - τ

Define equilibrium problem: set up equation system at each node and define iterates

State variable ω = Ω[[i]]:(c1/(c1+c2))=(c1/(δ[T-τ,j]))=Ω[[i]]    ⟺ c1=Ω[[i]]*δ[T-τ,j] &  c2=(1-Ω[[i]])*δ[T-τ,j]

Equations at each node

Clear[unknowns];unknowns[τ_,j_] := {c1u, c2u, c1d, c2d}; (*make sure to put braces even if there is only one unknown*)

equations[τ_, j_] := {(*Flow budget constraint or marketability condition*)
c2u + F2[τ - 1,j,c1u,c2u] == c2d + F2[τ - 1,j + 1,c1d,c2d], 
(*Kernel condition*)
((ϕ1[τ - 1, c1u] + ϕ1[τ - 1, c1d])/(ϕ1[τ, c1])) ==  ((ϕ2[τ - 1, c2u] + ϕ2[τ - 1, c2d])/(ϕ2[τ, c2])), 
   (*Aggregate resource constraints*)
 c1u + c2u == δ[T - τ + 1, j], c1d + c2d == δ[T - τ + 1, j + 1]};

(1 - Ω[[i]]) δ[T - τ, j]

initials[τ_, j_, i_] := {Ω[[i]]*δ[T - τ+1, j], (1 - Ω[[i]])*δ[T - τ+1, j], 
    Ω[[i]]*δ[T - τ + 1, j+1], (1 - Ω[[i]])*δ[T - τ + 1, j+1]};

Relate equation system to endogenous state variable

dummy[τ_, j_, i_] := {c1 -> Ω[[i]]*δ[T - τ, j], c2 -> (1 - Ω[[i]])*δ[T - τ, j]}

Complete - market equations

completeMarketEquations[τ_, j_] := {ϕ1[τ - 1, c1u]/ϕ1[τ, c1] == 
     ϕ2[τ - 1, c2u]/ϕ2[τ, c2], ϕ1[τ - 1, c1d]/ϕ1[τ, c1] == 
     ϕ2[τ - 1, c2d]/ϕ2[τ, c2], c1u + c2u == δ[T - τ + 1, j], 
    c1d + c2d == δ[T - τ + 1, j + 1]};

Define the function to be iterated backward, choosing a undiscountation for it that will facilitate the interpolation : below, state "undiscounted iterate" in the undiscounted form

F2[0,j_,c1_,c2_]:=0;
F2[τ_,j_,c1_,c2_]:=F2T[[τ,j,1]][(c1/(δ[T - τ, j]))]/;τ>0(*make sure to include a number in the third dimension of ϕ2F2T, even if there is only one iterate*)

UndiscountedIterate[τ_,j_,i_] := {((1/2)*(ϕ2[τ - 1, c2u]*(c2u+F2[τ - 1,j,c1u,c2u])) + (1/2)*(ϕ2[τ - 1, c2d]*(c2d+F2[τ - 1,j+1,c1d,c2d])))/ϕ2[τ,(1-Ω[[i]])δ[T - τ, j]]} ;(*make sure to have braces even if there is only one iterate*)

Input values

Number of dates and preference parameters

T = 6;   (* Note, first date = 0, so total number of dates is T + 1 *)
γ2=-5;
γ1:=-1;
ρ=1/(1+1/1000);

Total number of nodes at every point in time

TotK[t_]:=t+1;

Number of points for the center part of the grid: we use 2l+1 points

l = 50;

Precision goal for the endpoints

EndPointsGoal = 10^(-20);

Precision goal and working precision for FindRoot employed at each point

PG = 30; WP = 300;

The processes for output

δ0:= 100;
μδ=183/10000;σδ=357/10000;

Generate the process (binomial) for exogenous output

u := Exp[μδ-(((σδ)^2)/2)+σδ]-1;
d:=Exp[μδ-(((σδ)^2)/2)-σδ]-1;
δ[0,1] := δ0;
δ[t_, j_] := δ[t-1,Max[j-1,1]] *If[ j==1, (1 + u), (1+d)] ;

Size up the problem and introduce exponential form (one should not have to open this)

NbUnknowns[τ_,j_]:=Length[unknowns[τ,j]];

NbIterates:=Length[UndiscountedIterate[1,1,1] ]/.Table[unknowns[1,1]->initials[1, 1, 1] ];

Turn to exponential form to guarantee positivity of the solution

Exponentiate[τ_,j_] := Table[unknowns[τ, j][[ii]]->Exp[unknowns[τ, j][[ii]]],{ii,1,NbUnknowns[τ, j]}]

Eequations[τ_, j_] := equations[τ, j] /. Exponentiate[τ, j]

eUndiscountedIterate[τ_,j_,i_] :=UndiscountedIterate[τ,j,i]/. Exponentiate[τ, j];

Calculate function at time T - 1 (τ = 1), setting up ω - grid adaptively (one should not have to open this)

Center grid: we use 2l+1 points

F2Undiscount  will be a list of values of the Undiscounted function to be iterated backward at grid points. F2T will be a list of the same but interpolated.

F2T = {}; F2Undiscount = {}; Consumption = {};

τ = 1;Clear[j];

Ω = Table[(i - 1)/(2*l), {i, 2*l + 1}];

(prec = WP; locListAll = {}; UndiscountListAll = {}; PrevForecastLAll = {}; 
  PrevForecastRAll = {}; accLAll = {}; precLAll = {}; accRAll = {}; precRAll = {}; 
  n := Length[Ω]; m := Length[UndiscountList]; ll := Length[locList]; j = 1; 
Do[(*Do loop over all the nodes j of time T - 1*)
locList = {}; UndiscountList = {}; 
   (*Centerpoint for that node*)
 lastSol = FindRoot[equations[τ, j] /. dummy[τ, j, l + 1], 
      Transpose[{unknowns[τ, j], initials[τ, j, l + 1]}], PrecisionGoal -> PG, 
      WorkingPrecision -> prec,MaxIterations->20000];AppendTo[locList, unknowns[τ, j] /. lastSol]; 
    AppendTo[UndiscountList,  UndiscountedIterate[τ,j,l+1] /. lastSol]; ii = 1;
   Do[(*Do loop over all the grid points ii of one node*)
(*Prediction: Predict solution at the next grid point to the right, to serve as starting solution next time we call FindRoot: if three points are available, use an extrapolation. If not, use the initial values provided when defining the equation systsem*)
guess = If[ii < 3, initials[τ, j, l + 1 + ii], 
      Table[(Fit[{{Ω[[l + ii - 3]], locList[[2*ii - 4,k]]}, {Ω[[l + ii - 2]], 
               locList[[2*ii - 3,k]]}, {Ω[[l + ii - 1]], locList[[2*ii - 2,k]]}, 
              {Ω[[l + ii]], locList[[2*ii - 1,k]]}}, {1, y, y^2, y^3}, y] /. 
            {y -> Ω[[l + 1 + ii]]}), {k,NbUnknowns[τ, j]}]];
(*Correction: Calculate the next point to the right*)
      lastSol = FindRoot[equations[τ, j] /. dummy[τ, j, l + 1 + ii], 
        Transpose[{unknowns[τ, j], guess}], PrecisionGoal -> PG, WorkingPrecision -> prec,MaxIterations->20000];  (*Record the result: locList contains the results for equilibrium consumptions and UndiscountList contains the result for the Undiscounted function*)
 AppendTo[locList, unknowns[τ, j] /. lastSol]; AppendTo[UndiscountList, 
      UndiscountedIterate[τ,j,l+1+ii]/. 
        lastSol];
   (*Prediction: predict solution at the next grid point to the left, to serve as starting solution next time we call FindRoot: if three points are available, use an extrapolation. If not, use the initial values provided when defining the equation systsem*) 
guess = If[ii < 3, initials[τ, j, l + 1 - ii], 
Table[(Fit[{{Ω[[l + 2 - ii]], locList[[1,k]]}, {Ω[[l + 3 - ii]], locList[[
                2,k]]}, {Ω[[l + 4 - ii]], locList[[3,k]]}, {Ω[[l + 5 - ii]], 
               locList[[4,k]]}}, {1, y, y^2, y^3}, y] /. {y -> Ω[[l + 1 - ii]]}), 
           {k,NbUnknowns[τ, j]}]];
(*Correction: Calculate the next point to the right*)
lastSol = FindRoot[equations[τ, j] /. 
         dummy[τ, j, l + 1 - ii], Transpose[{unknowns[τ, j], guess}], PrecisionGoal -> PG, 
        WorkingPrecision -> prec]; 
  (*Record the result*)
PrependTo[locList, unknowns[τ, j] /. lastSol]; PrependTo[UndiscountList,   UndiscountedIterate[τ,j,l+1-ii]/. lastSol];
 ii = ii + 1; , {l - 1}]; 
   (*Record all the gridpoint values of the Undiscounted function for this node*) 
AppendTo[locListAll, locList]; AppendTo[UndiscountListAll, UndiscountList]; 
   (*Extrapolate to the endpoints to keep track of how much the endpoint changes as one expands the grid*)
(*Right endpoint*)PrevForecastR = Table[(Fit[{{Ω[[n - 4]], UndiscountList[[m - 3,iterate]]}, 
        {Ω[[n - 3]], UndiscountList[[m - 2,iterate]]}, {Ω[[n - 2]], UndiscountList[[m - 1,iterate]]}, 
        {Ω[[n - 1]], UndiscountList[[m,iterate]]}}, {1, y, y^2, y^3}, y] )/. {y -> Ω[[n]]},{iterate,1,NbIterates}];
   (*Record the resulting forecast into a list called "PrevForecastRAll"*)
 AppendTo[PrevForecastRAll, PrevForecastR];
(*Left endpoint*)PrevForecastL = 
    Table[(Fit[{{Ω[[2]], UndiscountList[[1,iterate]]}, {Ω[[3]], UndiscountList[[2,iterate]]}, 
        {Ω[[4]], UndiscountList[[3,iterate]]}, {Ω[[5]], UndiscountList[[4,iterate]]}}, {1, y, y^2, y^3}, 
       y] /. {y -> Ω[[1]]}),{iterate,1,NbIterates}];
  (*Record the resulting forecast into a list called "PrevForecastLAll"*)
AppendTo[PrevForecastLAll, PrevForecastL];
(*Initiate lists that will be used in the next subsection to keep track of the precision of the endpoints*) 
    AppendTo[accLAll, 1]; AppendTo[precLAll, 1]; AppendTo[accRAll, 1]; 
    AppendTo[precRAll, 1];
(*Go to the next node*)
 j = j + 1,{TotK[T-τ]}]; Clear[j]; )//Timing

For all nodes, expand Ω to the left until either the extrapolation precision at all nodes, or the extrapolation accuracy at all nodes becomes ⩽EndPointsGoal

As we get closer to the edges of the ω grid, we prefer to use the equation system in the exponential form, lest FindRoot attempt to compute with negative values of consumption*)

(While[Min[Max[precLAll], Max[accLAll]] > EndPointsGoal, 
   precLAll = {}; accLAll = {}; 
(*Replace the first point in the ω-list with one that is closer to the second point*)
Ω[[1]] = (1/2)*Ω[[2]]; PrependTo[Ω, 0]; prec = 1000; 
    j = 1; 
Do[locList = locListAll[[j]]; UndiscountList = UndiscountListAll[[j]]; 
      PrevForecastL = PrevForecastLAll[[j]];
      guess = Table[Fit[{{Ω[[3]], locList[[1,k]]}, {Ω[[4]], locList[[2,k]]}, 
           {Ω[[5]], locList[[3,k]]}, {Ω[[6]], locList[[4,k]]}}, {1, y, y^2, y^3}, y] /. 
         {y -> Ω[[2]]}, {k,NbUnknowns[τ, j]}]; lastSol = FindRoot[Eequations[τ, j] /. 
         dummy[τ, j, 2], Transpose[{unknowns[τ, j], Log[guess]}], PrecisionGoal -> PG, 
        WorkingPrecision -> prec, MaxIterations -> 2000]; PrependTo[locList, 
       Exp[unknowns[τ, j] /. lastSol]]; PrependTo[UndiscountList, 
      eUndiscountedIterate[τ,j,2]/. lastSol];LastForecastL = 
      Table[ Fit[{{Ω[[2]], UndiscountList[[1,iterate]]}, {Ω[[3]], UndiscountList[[2,iterate]]}, 
          {Ω[[4]], UndiscountList[[3,iterate]]}, {Ω[[5]], UndiscountList[[4,iterate]]}}, {1, y, y^2, y^3}, 
         y] /. {y -> Ω[[1]]},{iterate,1,NbIterates}];precL = Abs[Abs[LastForecastL]/Abs[PrevForecastL] - 1]; 
      accL = Abs[Abs[LastForecastL] - Abs[PrevForecastL]]; PrevForecastL = LastForecastL; 
      locListAll[[j]] = locList; UndiscountListAll[[j]] = UndiscountList; 
      PrevForecastLAll[[j]] = PrevForecastL; AppendTo[precLAll, precL]; 
      AppendTo[accLAll, accL]; j = j + 1; ,{TotK[T-τ]}]
]; Clear[j])//Timing

(*number of points added to the left*)n1 = n - 2*l - 1

For all nodes, expand Ω to the right until either the extrapolation precision at all nodes, or the extrapolation accuracy at all nodes becomes ⩽EndPointsGoal

(While[Min[Max[precRAll], Max[accRAll]] > EndPointsGoal, precRAll = {}; accRAll = {}; 
  (*Replace the last point in the ω-list with one that is closer to the last-but-one point*)
  Ω[[n]] = (1/2)*Ω[[n - 1]] + (1/2)*Ω[[n]]; AppendTo[Ω, 1]; prec = 1000; j = 1; 
    Do[locList = locListAll[[j]]; UndiscountList = UndiscountListAll[[j]]; 
      PrevForecastR = PrevForecastRAll[[j]]; 
      guess = Table[Fit[{{Ω[[n - 5]], locList[[ll - 3,k]]}, {Ω[[n - 4]], 
            locList[[ll - 2,k]]}, {Ω[[n - 3]], locList[[ll - 1,k]]}, 
           {Ω[[n - 2]], locList[[ll,k]]}}, {1, y, y^2, y^3}, y] /. {y -> Ω[[n - 1]]}, 
        {k,NbUnknowns[τ, j]}]; lastSol = FindRoot[Eequations[τ, j] /. dummy[τ, j, n - 1], 
        Transpose[{unknowns[τ, j], Log[guess]}], PrecisionGoal -> PG, WorkingPrecision -> prec, 
        MaxIterations -> 2000]; AppendTo[locList, Exp[unknowns[τ, j] /. lastSol]]; 
      AppendTo[UndiscountList, eUndiscountedIterate[τ,j,n-1] /. lastSol]; 
      LastForecastR = Table[Fit[{{Ω[[n - 4]], UndiscountList[[ll - 3,iterate]]}, 
          {Ω[[n - 3]], UndiscountList[[ll - 2,iterate]]}, {Ω[[n - 2]], UndiscountList[[ll - 1,iterate]]}, 
          {Ω[[n - 1]], UndiscountList[[ll,iterate]]}}, {1, y, y^2, y^3}, y] /. {y -> Ω[[n]]},{iterate,1,NbIterates}]; 
      precR = Abs[Abs[LastForecastR]/Abs[PrevForecastR] - 1]; 
      accR = Abs[Abs[LastForecastR] - Abs[PrevForecastR]]; PrevForecastR = LastForecastR; 
      locListAll[[j]] = locList; UndiscountListAll[[j]] = UndiscountList; 
      PrevForecastRAll[[j]] = PrevForecastR; AppendTo[precRAll, precR]; 
      AppendTo[accRAll, accR]; j = j + 1; , {TotK[T-τ]}]]; Clear[j])//Timing

Dimensions[ UndiscountListAll]

(*number of points added to the right*)n2 = n - n1 - 2*l - 1

N[{Max[precRAll], Max[accRAll]}]

The next piece of code does the following:
1) Checks if the final forecast for the value of F2T at omega=0 for node j is negative but can be treated as 0 within the precision specified by the variable EndPointsGoal;if this is the case then it sets the forecast to 0 – this ensures that the list of values that will be interpolated later will not contain 0;
2) does the same with regard to the forecast for of F2T at omega=1;
3) appends to the list “Consumption” the consumption levels calculated (for the present value of tau) at all nodes j (which are stored in the list “locListAll”);thus “Consumption[[tau]]” contains the consumption levels at all nodes at time tau,“Consumption[[tau,j]]” contains the consumption levels at time tau for node j,and so on.4) creates (for the current tau) the lists tmpList and tmpF2T so that tmpList[[j]] contains all the values that will be interpolated at node j (including the predicted and possibly corrected values at the left and the right endpoint) and tmpF2T[[j]] contains the actual interpolating function on the entire interval;
5) appends the lists tmpList and tmpF2T (which are local to period tau) to the global lists F2Undiscount and F2T;this has the effect that the Undiscounted function associated with period tau and node j is simply F2T[[tau,j]][omega] (tau and j are indices and only omega is an argument),while the values from which this function is constructed are stored in F2Undiscount[[tau,j]],which is itself a list.

(Do[If[0 > PrevForecastRAll[[j]] > EndPointsGoal, PrevForecastRAll[[j]] = 0], 
   {j, TotK[T-τ]}]; Do[If[0 > PrevForecastLAll[[j]] > EndPointsGoal, 
    PrevForecastLAll[[j]] = 0], {j,TotK[T-τ]}]; AppendTo[Consumption, locListAll]; 
  tmpList = {}; tmpF2T = {}; Do[loc = UndiscountListAll[[j]]; 
    PrependTo[loc, PrevForecastLAll[[j]]]; AppendTo[loc, PrevForecastRAll[[j]]]; 
    AppendTo[tmpList, loc]; AppendTo[tmpF2T, Table[Interpolation[Transpose[{Ω, loc[[All,iterate]]}]],{iterate,1,NbIterates}]], 
   {j,TotK[T-τ]}]; AppendTo[F2T, tmpF2T]; AppendTo[F2Undiscount, tmpList]; )//Timing

Do all remaining time iterations on the same grid (one should not have to open this)

τ=1;(Do[τ = τ + 1; prec = WP; locListAll = {}; UndiscountListAll = {}; PrevForecastLAll = {}; 
   PrevForecastRAll = {}; n := Length[Ω]; m := Length[UndiscountList]; 
   ll := Length[locList]; j = 1; 

Do[(*for each node*)locList = {}; UndiscountList = {}; 

     Do[(*for each point to the left of the center*)
guess = If[ii < n1 + l + 3, initials[τ, j, ii], 
         Table[Fit[{{Ω[[ii - 4]], locList[[ll - 3,k]]}, {Ω[[ii - 3]], 
                locList[[ll - 2,k]]}, {Ω[[ii - 2]], locList[[ll - 1,k]]}, {
                Ω[[ii - 1]], locList[[ll,k]]}}, {1, y, y^2, y^3}, y] /. 
             {y -> Ω[[ii]]}, {k, NbUnknowns[τ, j]}]]; 
       lastSol = FindRoot[equations[τ, j] /. dummy[τ, j, ii], 
         Transpose[{unknowns[τ, j], guess}], PrecisionGoal -> PG, WorkingPrecision -> prec]; 
       AppendTo[locList, unknowns[τ, j] /. lastSol]; AppendTo[UndiscountList, 
    UndiscountedIterate[τ,j,ii]/. lastSol]; , 
      {ii, n1 + l - 1, n - 1}
];
     Do[(*for each point to the right of the center*)
guess = ( Table[Fit[{{Ω[[ii + 1]], locList[[1,k]]}, {Ω[[ii + 2]], locList[[
                2,k]]}, {Ω[[ii + 3]], locList[[3,k]]}, {Ω[[ii + 4]], locList[[4,
                k]]}}, {1, y, y^2, y^3}, y] /. {y -> Ω[[ii]]}, {k, NbUnknowns[τ, j]}]); lastSol = FindRoot[equations[τ, j] /. dummy[τ, j, ii], 
         Transpose[{unknowns[τ, j], guess}], PrecisionGoal -> PG, WorkingPrecision -> prec, 
         MaxIterations -> 1000]; PrependTo[locList, unknowns[τ, j] /. lastSol]; 
       PrependTo[UndiscountList,UndiscountedIterate[τ,j,ii]/. lastSol]; , 
{ii, n1 + l - 2, 2, -1}
]; AppendTo[locListAll, locList]; AppendTo[UndiscountListAll, UndiscountList]; 
    
(*Add the righmost point*)
PrevForecastR = Table[
Fit[{{Ω[[n - 4]], UndiscountList[[m - 3,iterate]]}, 
         {Ω[[n - 3]], UndiscountList[[m - 2,iterate]]}, {Ω[[n - 2]], UndiscountList[[m - 1,iterate]]}, 
         {Ω[[n - 1]], UndiscountList[[m,iterate]]}}, {1, y, y^2, y^3}, y] /. {y -> Ω[[n]]},
{iterate,1,NbIterates}
];AppendTo[PrevForecastRAll, PrevForecastR]; 

(*Add the leftmost point*)
PrevForecastL = 
      Table[
Fit[{{Ω[[2]], UndiscountList[[1,iterate]]}, {Ω[[3]], UndiscountList[[2,iterate]]}, 
         {Ω[[4]], UndiscountList[[3,iterate]]}, {Ω[[5]], UndiscountList[[4,iterate]]}}, {1, y, y^2, y^3}, 
        y] /. {y -> Ω[[1]]},
{iterate,1,NbIterates}
]; AppendTo[PrevForecastLAll, PrevForecastL]; 

j = j + 1; , 
    {TotK[ T-τ]}
];

(*For each node, replace extremeties with negative values*)
 Clear[j]; Do[
If[0 > PrevForecastRAll[[j]] > EndPointsGoal, 
     PrevForecastRAll[[j]] = 0], {j, TotK[ T-τ]}]; Do[
If[0 > PrevForecastLAll[[j]] > EndPointsGoal, PrevForecastLAll[[j]] = 0], 
    {j, TotK[ T-τ]}
]; AppendTo[Consumption, locListAll]; 

(*For each node, interpolate*)
tmpList = {}; tmpF2T = {}; 
   Do[
loc = UndiscountListAll[[j]]; PrependTo[loc, PrevForecastLAll[[j]]]; 
     AppendTo[loc, PrevForecastRAll[[j]]]; AppendTo[tmpList, loc]; 
     AppendTo[tmpF2T, Table[Interpolation[Transpose[{Ω, loc[[All,iterate]]}]],{iterate,1,NbIterates}]], 
{j, TotK[ T-τ]}
]; AppendTo[F2T, tmpF2T]; AppendTo[F2Undiscount, tmpList]; 

, 
{T-1}
])//Timing

Generate output at time 0

τ = T;

SecondPlot=Plot[Table[F2T[[τ,j,1]][1-ω], {j, T - τ + 1}], {ω, 0, 1},PlotRange->All,Frame->True,FrameTicks->Automatic,FrameTicksStyle->Directive[Black,10],PlotLabel->Style[Column[{"Non Stockholders' Wealth at time 0"}],FontSize->15],FrameLabel->{Style["Non Stockholders' Share of Aggregate Consumption at time 0",FontSize->15],},RotateLabel->False,PlotStyle->{AbsoluteThickness[3]}]

τ = T; j = 1;FirstPlot= ListPlot[
   {Table[{1 - Ω[[i + 1]], (1/2)*(ϕ1[τ - 1, ConsumptionCM[[i,3]]]/
         ϕ1[τ, Ω[[i + 1]]*δ[T - τ, j]] - ϕ1[τ - 1, ConsumptionCM[[i,1]]]/
         ϕ1[τ, Ω[[i + 1]]*δ[T - τ, j]])}, {i, m}], 
    Table[{1 - Ω[[i + 1]], (1/2)*(ϕ1[τ - 1, Consumption[[τ,j,i,3]]]/
         ϕ1[τ, Ω[[i + 1]]*δ[T - τ, j]] - ϕ1[τ - 1, Consumption[[τ,j,i,1]]]/
         ϕ1[τ, Ω[[i + 1]]*δ[T - τ, j]])}, {i, m}]}, PlotRange->All,Frame->True,FrameTicks->Automatic,FrameTicksStyle->Directive[Black,10],PlotLabel->Style[Column[{"Market Price of Risk"}],FontSize->15],Frame->True,FrameTicks->Automatic,FrameTicksStyle->Directive[Black,20],FrameLabel->{Style["Non Stockholders' Share of Aggregate Consumption at time 0",FontSize->15]},RotateLabel->False,PlotMarkers->{{"■",3},{"▼",8}},Epilog->Inset[Style[Column[{"Incomplete ▼ ▼ ▼ ▼ ▼","Complete    ◼ ◼ ◼ ◼ ◼"}],15],{0.2,0.8},Automatic,{0.4,0.3}]]

Show[GraphicsColumn[{FirstPlot,SecondPlot}]]

Export["D:\Mes documents\Rech\Incomple\MktPriceofRisk_not_reverse.eps",%, ImageSize-> {72*10, Automatic}]

Show[GraphicsRow[{FirstPlot,SecondPlot}]]

Export["D:\Mes documents\Rech\Incomple\MktPriceofRisk_not_reverse_row.eps",%, ImageSize-> {72*10, Automatic}]

τ=T;j=1;ii=75;ContourPlot[Evaluate[(equations[τ, j][[1;;2]] )/.{c1u->-c2u+δ[T - τ + 1, j],c1d->-c2d+δ[T - τ + 1, j+1]}/. dummy[τ, j, ii]],{c2u,1/100,δ[T - τ + 1, j]-1/100},{c2d,1/100,δ[T - τ + 1, j+1]-1/100},WorkingPrecision->100,FrameLabel->{"c2u","c2d"},LabelStyle->Directive[Black,FontSize->20],ContourLabels->({Text[Style["Kernel condition",FontSize->25],{70,27},Automatic,{1,-1/3}],Text[Style["Marketability condition",FontSize->25],{70,75},Automatic,{1,1}],Text[Style["Equilibrium",FontSize->25],{30,42}]}&),ContourStyle->{AbsoluteThickness[2]}]

Export["D:\Mes documents\Rech\Incomple\geometry.eps",%, ImageSize-> {72*10, Automatic}]