This file implements the algorithm 'Forward Shoot 2D Anticipated'  outlined in paper titled, "Smart Forward Shooting" which solves by forward shooting the models with two jump variables when the shock is anticipated.

The algorithm and the program written to implement the algorithm is free for use as long as proper credit is given to the authors of the paper and this first cell of the program is kept unchanged. This program is licensed under the CC-GNU GPL license with no expressed or implied guarantees.

Works with Mathematica Version 5.2 or later.

Version: April 07, 2008.
Author: Manoj Atolia

For details of the particular model being solved here, see Example 5 of the paper.

Off[General::spell];Off[General::spell1];Off[Unset::norep];Off[Unset::wrsym];Off[InterpolatingFunction::dmval];Off[Precision::precsm];

(* The precision of preliminary calculations. *)
{precision = 100, $MinPrecision=digits=precision};

(* CHANGES WITH OFFICIAL VERSION:
  1. REMOVED SESSION TIME COMMANDS IN OFFICIAL VERSION. 
2. CHANGED THE BOX SIZE TO NORMAL IN OFFICIAL VERSION. THE BOXES ARE MUCH BIGGER IN ANTICIPATED CASE.
3. THERE IS NO SHRINKING OF RING ON FLIPPING IN THE OFFICIAL VERSION. SEE THE RELEVANT COMMAND COMMENTED OUT IN COMDIS2.	*)
(* FURTHER IMPROVEMENTS IN FUTURE:
  1. DYNAMICALLY CHANGING THE BOX. INITIALLY STARTING WITH POSITIVE PART OF SPACE AND THEN MOVING TO A REASONABLE BOX. MORE RELEVANT FOR TEMPORARY AND ANTICIPATED CASES AS WELL AS TO HIGHER DIMENSIONS.
2. INTERACTIVE INPUT WITH Input[] COMMAND. *)

Initialization

(* Set the directory where Mathematica will store intermediate results. *)
SetDirectory["c:"];

(* ALWAYS ENTER NUMERICAL VALUES AS FRACTIONS AND NOT AS DECIMALS TO ENSURE REQUISITE PRECISION IN CALCULATIONS *)   
{τ=(1/2),σ=(1/3),ρ=(1/10),β=(1/2),δ=(1/20)}//N                                   (* Parameters of the model *)

{(Subscript[bbyY, o])=(25/100),(Subscript[mbyY, o])=(8/100),(Subscript[π, o])=(25/100),(Subscript[θ, o])=(1/2)}//N        (* Data for calibration *)

{(Subscript[k, o])=1}//N                                                                          (* Normalization *)

{(Subscript[v, f])=(15/100),(Subscript[v, o])=(20/100),(Subscript[ω, 1])=(1/4),(Subscript[ω, 2])=(1/2),ts=2}//N                                                                 (* Shock - Exogenous rise in productivity *)

Steady State

(* Define utility function and production function *)
U[c_,m_]:=((1 )/(1-(1/τ)))(((((((c^(1-(1/σ)))+g (m^(1-(1/σ)))))^(σ/(σ-1)))))^(1-(1/τ)));(Subscript[U, 1])[c_,m_]:=Evaluate[D[U[c,m],c]];(Subscript[U, 2])[c_,m_]:=Evaluate[D[U[c,m],m]];
Y[k_]:=((((Subscript[a, 0])+(Subscript[a, 1])(k^((β-1)/β))))^(β/(β-1)));(Subscript[Y, 1])[k_?NumberQ]:=Evaluate[D[Y[k],k]];
(* Avoid using symbols for subscripts for representing derivatives. Thus, (Subscript[Y, 1])[k,z] is (Subscript[Y, K])[k,z] of the paper.*)
(* Since we are using symbols o and f as subscripts for initial and final steady states avoid using o and f as names of variables or function. *)

(* Calibrate and compute the original steady state *)
(Subscript[Y, o])=Y[(Subscript[k, o])];{(Subscript[c, o])=(Subscript[Y, o])-δ (Subscript[k, o]),(Subscript[m, o])= (Subscript[mbyY, o]) (Subscript[Y, o]),b=(Subscript[bbyY, o]) (Subscript[Y, o]),(Subscript[z, o])=(Subscript[π, o])(Subscript[m, o])+(Subscript[v, o])Y[(Subscript[k, o])]-ρ b};g=((Subscript[π, o])+ρ)(((((Subscript[m, o])/(Subscript[c, o]))))^(1/σ));{(Subscript[a, 0]),(Subscript[a, 1])}=.;{(Subscript[a, 0]),(Subscript[a, 1])}={(Subscript[a, 0]),(Subscript[a, 1])}/.FindRoot[{(Subscript[Y, 1])[(Subscript[k, o])] == ((ρ+δ)/(1-(Subscript[v, o]))),(Subscript[θ, o])==(((ρ+δ)(Subscript[k, o]))/(Subscript[Y, o]))},{(Subscript[a, 0]),1},{(Subscript[a, 1]),1},MaxIterations->1000,AccuracyGoal->digits,WorkingPrecision->digits];

(* Compute the new steady state *)
(Subscript[k, f])=.;(Subscript[k, f])=(Subscript[k, f])/.FindRoot[{(Subscript[Y, 1])[(Subscript[k, f])]==((ρ+δ)/(1-(Subscript[v, f])))},{(Subscript[k, f]),(Subscript[k, o])},MaxIterations->1000,AccuracyGoal->digits,WorkingPrecision->digits];(Subscript[Y, f])=Y[(Subscript[k, f])];(Subscript[c, f])=(Subscript[Y, f])-δ (Subscript[k, f]);(Subscript[z, f])=(Subscript[z, o])+(Subscript[v, f])Y[(Subscript[k, f])]-(Subscript[v, o])Y[(Subscript[k, o])] ;
{(Subscript[m, f]),(Subscript[π, f])}=.;{(Subscript[m, f]),(Subscript[π, f])}={(Subscript[m, f]),(Subscript[π, f])}/.FindRoot[{g(((((Subscript[m, f])/(Subscript[c, f]))))^(-1/σ))==(Subscript[π, f])+ρ,(Subscript[z, f])==(Subscript[π, f])(Subscript[m, f])+(Subscript[v, f])Y[(Subscript[k, f])]-ρ b},{(Subscript[m, f]),(3/2)(Subscript[m, o])},{(Subscript[π, f]),(Subscript[π, o])},MaxIterations->1000,AccuracyGoal->digits,WorkingPrecision->digits];
{{m,c,k,Y,"b",z,"pi","c/m","m/Y"},{(Subscript[m, o]),(Subscript[c, o]),(Subscript[k, o]),(Subscript[Y, o]),b,(Subscript[z, o]),(Subscript[π, o]),(Subscript[c, o])/(Subscript[m, o]),(Subscript[m, o])/(Subscript[Y, o])},{(Subscript[m, f]),(Subscript[c, f]),(Subscript[k, f]),(Subscript[Y, f]),b,(Subscript[z, f]),(Subscript[π, f]),(Subscript[c, f])/(Subscript[m, f]),(Subscript[m, f])/(Subscript[Y, f])}}//N//MatrixForm

Global Nonlinear Saddle Path

(* Specify the nonlinear system with jump variables in the core dynamic system in the beginning and the state variables at the end *)
{
mdot = z+r b-v Y[k] - pi m; 
kdot = Y[k]-c-δ k;
λdot =(ρ-((1-v)(Subscript[Y, 1])[k]-δ))λ;
cdot=dcin/.Solve[Dt[(Subscript[U, 1])[c,m]==λ]/.{Dt[c]->dcin,Dt[m]->mdot,Dt[k]->kdot,Dt[v]->vdot,Dt[λ]->λdot},dcin][[1]]/.λ->(Subscript[U, 1])[c,m];
system  ={mdot,cdot,kdot};
};

(* Solve for the endogenous variables as functions of the variables in the system. For consistency maintain same order of variables in arguments for functions of endogenous variables as in the system *)
{
rfunc[k_,v_]:=(1-v)(Subscript[Y, 1])[k]-δ;
pifunc[m_,c_, k_,v_]:=g((((m/c)))^(-(1/σ)))-rfunc[k,v];
};

(* Specify the functions, IF ANY, that, for a given value of t, return the values of the (shock) variables that change exogenously. *)
z[t_]:=(Subscript[z, o])+((Subscript[z, f])-(Subscript[z, o]))(1-Exp[-(Subscript[ω, 1])(t-ts)])UnitStep[t-ts];
v[t_]:=(Subscript[v, o])+((Subscript[v, f])-(Subscript[v, o]))(1-Exp[-(Subscript[ω, 2])(t-ts)])UnitStep[t-ts];

(* Specify the function that returns a list containing the "post-shock" values (the values of all variables in the system just after the impact of the anticipated shock) when given their "pre-shock" values (the values just before the impact of the shock).*)
func[m_,c_,k_,v_]:={m,c,k,v};

(* INITIALIZATION CELL - Initialization for solving for the global saddle path. *)
{  
T = 20;                        (* solution time (years). BEWARE!!! Setting T larger than necessary can substantially increase computational time.*); 
tol = (10^(-4));                (* desired accuracy of the solutions of the variables in the system *);
precision = 16;      (* starting value of precision of calculations*) ;
epsilon= ((10^(-0))/40)(((((Subscript[m, o])^2)+((Subscript[c, o])^2)))^(1/2));           (* radius of the ϵ-ring*)
pini={(Subscript[m, o]),(Subscript[c, o])};                    (* center of the initial ϵ-ring*);
};
(* IN THE UNLIKELY EVENT THAT THE PROGRAM FAILS AT THE FIRST ATTEMPT TRY INCREASING THE VALUE OF 'precision'. ONE MAY ALSO EXPERIMENT WITH DIFFERENT VALUES OF 'epsilon' and 'pini'. *)

(* EXECUTE THIS CELL TO SEE THE SYNTAX FOR THE solve COMMAND *)
?solve

(* MAIN CELL - Solving for the global nonlinear saddlepath. *)
{gtmin,gminsol1,gminsol2}=solve[system,{t},{m,c,k,v},{(Subscript[m, o]),(Subscript[c, o]),(Subscript[k, o]),(Subscript[v, o])},{(Subscript[m, f]),(Subscript[c, f]),(Subscript[k, f]),(Subscript[v, f])},{{r,rfunc,{k,v}},{pi,pifunc,{m,c,k,v}},{z,z,{t}}},{},{},{ts,func},{T,tol,precision,epsilon},pini]
Print["Maximum amount of memory used in computing the solution: ",MaxMemoryUsed[]/1000000MB//N]

(* Plot the final solutions. *)
tplot = T+5;
pc=Plot[{If[t<0,(Subscript[c, o]),c[t]]/.If[t<ts,gminsol1,gminsol2],(Subscript[c, f])},{t,-tplot/10,tplot},PlotRange->{.245,.262},AxesOrigin->{-tplot/10,.245},AxesLabel->{"Years",None},PlotStyle->{{Thickness[.007],Dashing[{.05,.0}]},{Thickness[.007],Dashing[{.02,0.02}]},{Thickness[.007],Dashing[{.01,0.01}]}},PlotLabel->"Real Consumption",DisplayFunction->Identity];
pm=Plot[{If[t<0,(Subscript[m, o]),m[t]]/.If[t<ts,gminsol1,gminsol2],(Subscript[m, f])},{t,-tplot/10,tplot},PlotRange->{.0188,.026},AxesOrigin->{-tplot/10,.0188},AxesLabel->{"Years",None},PlotStyle->{{Thickness[.007],Dashing[{.05,.0}]},{Thickness[.007],Dashing[{.02,0.02}]},{Thickness[.007],Dashing[{.01,0.01}]}},PlotLabel->"Real Money Balances",DisplayFunction->Identity];
pk=Plot[{If[t<0,(Subscript[k, o]),k[t]]/.If[t<ts,gminsol1,gminsol2],(Subscript[k, f])},{t,-tplot/10,tplot},PlotRange->{.99,1.09},AxesOrigin->{-tplot/10,.99},AxesLabel->{"Years",None},PlotStyle->{{Thickness[.007],Dashing[{.05,.0}]},{Thickness[.007],Dashing[{.02,0.02}]},{Thickness[.007],Dashing[{.01,0.01}]}},PlotLabel->"Capital Stock",DisplayFunction->Identity];
pY=Plot[{If[t<0,Y[(Subscript[k, o])],Y[k[t]]]/.If[t<ts,gminsol1,gminsol2],Y[(Subscript[k, f])]},{t,-tplot/10,tplot},PlotRange->{.298,.316},AxesOrigin->{-tplot/10,.298},AxesLabel->{"Years",None},PlotStyle->{{Thickness[.007],Dashing[{.05,.0}]},{Thickness[.007],Dashing[{.02,0.02}]},{Thickness[.007],Dashing[{.01,0.01}]}},PlotLabel->"Real Output",DisplayFunction->Identity];
pr=Plot[{If[t<0,rfunc[(Subscript[k, o]),(Subscript[v, o])],rfunc[k[t],v[t]]]/.If[t<ts,gminsol1,gminsol2],r[(Subscript[k, f]),(Subscript[v, f])]},{t,-tplot/10,tplot},PlotRange->{.098,.106},AxesOrigin->{-tplot/10,.098},AxesLabel->{"Years",None},PlotStyle->{{Thickness[.007],Dashing[{.05,.0}]},{Thickness[.007],Dashing[{.02,0.02}]},{Thickness[.007],Dashing[{.01,0.01}]}},PlotLabel->"Real Interest Rate",DisplayFunction->Identity];
ppi=Plot[{If[t<0,pifunc[(Subscript[m, o]),(Subscript[c, o]),(Subscript[k, o]),(Subscript[v, o])],pifunc[m[t],c[t],k[t],v[t]]]/.If[t<ts,gminsol1,gminsol2],pi[(Subscript[m, o]),(Subscript[c, o]),(Subscript[k, o]),(Subscript[v, f])]},{t,-tplot/10,tplot},PlotRange->{.20,.57},AxesOrigin->{-tplot/10,.20},AxesLabel->{"Years",None},PlotStyle->{{Thickness[.007],Dashing[{.05,0}]},{Thickness[.007],Dashing[{.02,0.02}]},{Thickness[.007],Dashing[{.01,0.01}]}},PlotLabel->"Inflation",DisplayFunction->Identity];

If[$VersionNumber==5.2,
Show[GraphicsArray[{{pc,pm},{pY,pk},{pr,ppi}}]],
Show[GraphicsGrid[{{pc,pm},{pY,pk},{pr,ppi}}]]]

(* Tabulate the final solutions. *)
Labels={"m","c","i","k","z","v","Y","r","pi"};InitialSteadyState={(Subscript[m, o]),(Subscript[c, o]),δ (Subscript[k, o]),(Subscript[k, o]),(Subscript[z, o]),(Subscript[v, o]),(Subscript[Y, o]),ρ,(Subscript[π, o])}//N;FinalSteadyState={(Subscript[m, o]),(Subscript[c, o]),δ (Subscript[k, o]),(Subscript[k, o]),(Subscript[z, o]),(Subscript[v, o]),(Subscript[Y, o]),ρ,(Subscript[π, o])}//N;
(Dat=Table[{tt,m[t],c[t],Y[k[t]]-c[t],k[t],z[t],v[t],Y[k[t]],rfunc[k[t],v[t]],pifunc[m[t],c[t],k[t],v[t]]}/.tt->t/.If[t<ts,gminsol1,gminsol2],{t,0,T,1/5}])//MatrixForm;
(Dat=Join[{Join[{"Initial Steady State"},InitialSteadyState],Join[{"Final Steady State"},FinalSteadyState],Join[{Time},Labels]},Flatten[Dat,1]])//N//MatrixForm

The solve forward shoot 2D Anticipated Shock Package

(* Interface *)
BeginPackage["solve`"];

EndPackage[];

Off[General::spell];Off[General::spell1];Off[Unset::norep];Off[Unset::wrsym];Off[InterpolatingFunction::dmval];Off[Precision::precsm];

(* Usage messages for the exported functions and the context itself *)
solve::usage = "solve[sys, {t}, vars, varso, varsf, varsendo, sysexpand, varspecial, antishock, init, pini] solves a system of first-order differential equations subjected to an anticipated shock at time ts > 0. \n
sys   = {xdot, ydot,...} lists the right hand side of the system of first-order differential equations (with jump variables in the core dynamic system in beginning and state variables in the end).\n
vars  = {x, y, ...,z} is the list of variables in the system in the same order as the differential equations in sys and the variable z changing in anticipated manner appended at the end. Note that although z appears in the list there is not differential equation corresponding to it. solve assumes that the values of z for different t are given by a function of same name whose only argument is t.\n
varso = {(Subscript[x, o]), (Subscript[y, o]), ...,(Subscript[z, o])} is the list of original steady values of variables in the system in the same order as vars.\n
varsf = {(Subscript[x, f]), (Subscript[y, f]), ...,(Subscript[z, f])} is the list of final steady values of variables in the system in the same order as vars.\n
varsendo = {{w, wfunc, {x,y, ...}}, ...} is the list of functions {wfunc, ...} for the endogenous variables {w, ...} where wfunc is function of {x,y, ...} and {x,y,...} are in the order in which they appear as arguments of wfunc. If no endogenous variables in the system, enter {{,{}}}.\n
sysexpand  = {funcj, {x, y, ...}} contains the name of function funcj that takes values of variables in core dynamic
	system {x, y, ...} at time 0 as argument and returns a list containing values of other jump variables (in the same order as they appear in vars) that have been added to core dynamic system to obtain the expanded system sys. If core dynamic system has not been expanded, enter {}.\n
varspecial  = {v,...} is the list of variables that can potentially take non-positive values. It is assumed that other variables in vars have been defined so that they only take positive values. If no variable can potentially take non-positive values, enter {}.\n
antishock = {ts, func} gives the time ts at which anticipated shock hits the economy and func is the function with {x,y,...,z} as arguments that returns a list containing post-shock values of all variables in vars at time ts (in same order as the variables are listed in vars) when given their pre-shock values at ts-.\n
init  = {T, tol, precision, epsilon} sets the required duration of the solution to T, the required accuracy of solution to tol, starting precision to precision, and radius of initial ϵ-ring to epsilon. Note that setting T larger than necessary can significantly increase computation time. \n
pini  = {a, b} gives center of the initial ring in the jump space. The typical value is {(Subscript[x, o]), (Subscript[y, o])}, the initial steady state values of the jump variables in the core dynamic system. A better guess may somewhat help reduce computational time but is not necessary.";

(* Error messages for the exported objects *)
{  solve::sysDESys = "You called `1` with first argument `2` (the system of differential equations) that has variables other than those defined in `3` and `4`.";

solve::sysLength = "You called `1` with lengths of `2`, `3`, `4` and `5` not equal.";

solve::vars = "You called `1` with third argument `2` that is not a list of variables.";

solve::origVals = "You called `1` with fourth argument `2` that is not a list of numbers.";

solve::finVals = "You called `1` with fifth argument `2` that is not a list of numbers.";

solve::prec = "You called `1` with fourth and fifth arguments `2` and `3` that have numbers with precision less than 100.";

solve::endoVars = "You called `1` with sixth argument `2` that is not a list of lists of the form {w, wfunc, {x,y, ...}} where w is an endogenous variable and {x,y, ...} are the variables in the system as defined in `3`.";

solve::initVals = "You called `1` with last argument `2` that is not a list of numbers.";

solve::sysError = "You called `1` with a system that does not have three negative eigenvalues when linearized."};

(* Implementation *)
Begin["solve`Private`"];

(* Definition of the exported functions *)
solve[asys_List,{at_},avars_List,avarso_List,avarsf_List,avarsendo_List,asysexpand_List,avarspecial_List,aantishock_List,ainit_List,apini_List]:=Block[{sys,t,vars,varso,varsf,varsendo,init,antishock,sysexpand,varspecial,T,tol,precision,ts,Tmax,firstsol,precisionini,system,boxvars,box,A,localmin,pmin,epsilon,digits,ntrajini,gmindis,gtmin,gminsol1,gminsol2,initsol,ini,dis,tsol,tdis,doneflag},

{sys,{t},vars,varso,varsf,varsendo,sysexpand,varspecial,antishock,init,pmin}={asys,{at},avars,avarso,avarsf,avarsendo,asysexpand,avarspecial,aantishock,ainit,apini};
{T,tol,precision,epsilon}=init;ts=First[antishock];Tmax=4000;firstsol = 1;precisionini=precision;

(* Check input *)
checkInput[];

(* Make the system from input data and create box over which to solve the system *)
system=makeSys[];
boxvars=Map[Part[#,Flatten[Position[Map[MemberQ[varspecial,#]&,Drop[vars,-1]],False]]]&,{Drop[vars,-1],Drop[varso,-1],Drop[varsf,-1]}];
box:=Evaluate[Apply[Or,Join[MapThread[First[Map[(#1),{t}]]<((Min[#2,#3])/20)&,boxvars],MapThread[First[Map[(#1),{t}]]>10Max[#2,#3]&,boxvars]]]];

(* SOLVE FOR THE GLOBAL NONLINEAR SADDLEPATH *);
While[True,
localmin=0;

Label[begin];
$MinPrecision=digits=precision;ntrajini=20;gmindis=tol;
If[firstsol==0,{pmin,epsilon}=(<<solution2)[[{5,6}]]];

(* Calculate the minimum distance for the trajectory of the nonlinear system that starts from the center of the intial ϵ-ring. *)
{gmindis,gtmin,gminsol1,gminsol2,pmin}=compdis[pmin,0];

Print["AT THE CENTER OF THE INITIAL RING: ","    Minimum Distance from Final Steady State = ",gmindis//N,"    Time of Closest Approach = ",gtmin//N];

(* SOLVE - Solve for the global nonlinear saddle path *);
If[localmin==1,localmin=0;pmin=pmin+epsilon/10{1/2,1/2}];

{gminsol1,gminsol2,gtmin,gmindis,pmin,epsilon}=compdis2[pmin,epsilon];

If[gtmin==(10^100),
Print["Either the system is highly ill-conditioned or you might have encountered a local minima. Perturing the center of the initial ring and resetting 'precision' to its initial value"];precision=precisionini;localmin=1;Goto[begin]];
If[doneflag==1,Goto[end]];

];
 
(* DURATION CHECK - Check the duration for which the solution is accurate within tolerance 'tol'.*);
Label[end];
durationCheck[];

{gtmin,gminsol1,gminsol2}];

checkInput[___]:=Block[{},

If[Not[Length[sys]==(Length[vars]-1)==(Length[varso]-1)==(Length[varsf]-1)],
	Message[solve::sysLength,"'solve'","'sys'","'vars'","'varso'","'varsf'"];Abort[]];

If[Not[Length[vars]==Length[Cases[vars,_Symbol]]],Message[solve::vars,"'solve'","'vars'"];Abort[]];

If[Not[Length[varso]==Length[Cases[varso,_?NumericQ]]],Message[solve::origVals,"'solve'","'varso'"];Abort[]];

If[Not[Length[varsf]==Length[Cases[varsf,_?NumericQ]]],Message[solve::finVals,"'solve'","'varsf'"];Abort[]];

If[Min[Map[Precision[#]&,Join[varso,varsf]]]<100,Message[solve::prec,"'solve'","'varso'","'varsf'"];Abort[]];

If[Not[Length[varsendo]==Length[Cases[varsendo,{_Symbol,_Symbol,_List}]]]||Not[Length[Flatten[varsendo]]==Length[Cases[Flatten[varsendo],_Symbol]]],
	Message[solve::endoVars,"'solve'","'varsendo'","'vars'"];Abort[];];

If[Not[Length[init]==Length[Cases[init,_?NumericQ]]],Message[solve::initVals,"'solve'","'init'"];Abort[]];

];

makeSys[___]:=Block[{system},

system=sys/. Thread[MapThread[(#1)&,Transpose[varsendo]]->
					MapThread[ToExpression[StringJoin[ToString[ #2],"[",ToString[#3],"]"]] /.List->Sequence&, Transpose[varsendo]],
						List];

If[Length[Cases[system/.t->Infinity/.Thread[vars->varsf,List],_?NumberQ]]≠(Length[vars]-1),
	Message[solve::sysDESys,"'solve'","'sys'","'vars'","'varsendo'"];Abort[];];

system=Thread[Map[ToExpression[StringJoin[ToString[#],"'","[",ToString[t],"]"]]&,Drop[vars,-1]]==  
			(system=system/.Thread[vars->Map[ToExpression[StringJoin[ToString[#],"[",ToString[t],"]"]]&,vars],List]),
				List];

system];

distance[sol1_,sol2_,tsol_,tini_]:=Block[{dis,dislast,tdis,steps,tbeg,tend,intsz,traj,traj1,traj2,norm,norm2,temp},
dis= (10^7);dislast=(10^8);If[firstit==1,steps=Max[4Floor[tsol+1],450] ,steps=Max[4Floor[tsol+1-tini],50] ];
tbeg=tini;tend=tsol;intsz=tend-tbeg;temp=Last[varsf];
norm=Table[vars[[i]][t]-varsf[[i]],{i,1,Length[vars]-1}];
norm2=ToExpression[StringJoin[ToString[Last[vars]],"[",ToString[t],"]"]]-temp;

While[ (dislast-dis)>gmindis (10^(-7)),
	dislast=dis;
	traj1=Table[(norm/.t->s/.If[s<ts,sol1,sol2])[[1]],{s,tbeg,tend,(intsz/steps)}];
	traj2=1000Table[({If[s<ts ,norm2/.t->s/.sol1,norm2/.t->s/.sol2]})[[1]],{s,tbeg,tend,(intsz/steps)}];
	traj = Map[Norm[#,2]&,MapThread[Join[#1,#2]&,{traj1,traj2}]];
	dis=Min[traj];tdis=tbeg+(Position[traj,dis][[1,1]]-1)(intsz/steps);
	tbeg=Max[tdis-(intsz/steps),tbeg];tend=Min[tdis+(intsz/steps),tsol];intsz=tend-tbeg;steps=40 ;
];

{dis,tdis}];

(* Set up for solving for the global nonlinear saddle path *)
compdis2[p_,ϵ_]:=Block[{rmin,rtmin,rminsol1,rminsol2,rpmin,difftmin,notsatistmin,satistmin,diffprec,notsatisprec,satisprec,increaseprec,congratulations,done,step2,step3,rit,iflag,shrink,pmin,epsilon},

iflag=0;trajtrack=0;rit=1;shrink=False;{pmin,epsilon}={p,ϵ};

(* Messages *)
difftmin:=Print["Difference in values of variables between two paths with (Subscript[t, min])s differing by at least 5 years = ",sysDiff[0,T]//N];
notsatistmin="The (Subscript[t, min]) is NOT large enough. The variables changed by greater than 'tol' due to latest increase in (Subscript[t, min]).";
satistmin= "THE ACCURACY CHECK PASSED SUCCESSFULLY. The variables changed by less than 'tol' due to latest increase in (Subscript[t, min]).";
diffprec:=Print["Difference in values of variables due to latest increase in precision = ", sysDiff[0,T]//N];
notsatisprec="The precision is NOT satisfactory. The variables changed by greater than 'tol' due to latest increase in 'precision'.";
satisprec="THE PRECISION IS SATISFACTORY AS WELL. The variables changed by less than 'tol' due to latest increase in 'precision'.";increaseprec:=Print["Increased 'precision' by 4 digits to ",precision," digits."];

step2="Starting Step 2: Checking for the accuracy of the candidate solution.";
step3="Starting Step 3: Checking for the precision of the candidate solution.";

congratulations= "CONGRATULATIONS! The variables changed by less than 'tol' due to latest increase in 'precision'.";
done = "You are done! Hope to see you in another model. Continue with the next cell.";

If[firstsol==0,rit=0;Goto[errorcheck]];
Label[begin];
Print[""];Print["BEGINNING SEARCH ON RING ",rit," WITH CENTER ",SetPrecision[pmin,precision]];

{rmin,rtmin,rminsol1,rminsol2,rpmin}=compdis[pmin,epsilon];

If[rtmin==(10^100)&& precision≤80,
	precision=precision+4;$MinPrecision=digits=precision;Print["Increased 'precision' by 4 digits to ",precision," digits."];Goto[begin],
	If[rtmin==(10^100),gtmin=rtmin;Goto[end] ]];

If[rmin≤gmindis,
	gmindis=rmin;gtmin=rtmin;gminsol1=rminsol1;gminsol2=rminsol2;pmin=rpmin;
	If[shrink==True,epsilon=epsilon/2,epsilon=epsilon/1],
	(*If[(rit≠1&&iflag==1)||gtmin>ts+5||firstsol==0,shrink=True];*)(* NO SHRINKING ON FLIPPING IN ANTICIPATED CASE. *)
	epsilon=epsilon/2];

Print["AT THE END OF SEARCH ON RING ",rit,":    CURRENT MINIMUM DISTANCE = ",gmindis//N,"     TIME OF CLOSEST APPROACH = ",gtmin];

(* ERROR CHECK - Check for the accuracy of the solution over the time period of interest. *) 
Label[errorcheck];
If[firstsol==1,
	If[gtmin>T+5 && iflag==0,iflag=1;Print[SessionTime[]];{gminsol1,gminsol2,gtmin,gmindis,pmin,5epsilon}>>solution1;Print[""];Print[step2];];
	If[gtmin>Max[T+5 ,(<<solution1)[[3]]+5]&& iflag==1 ,
		{gminsol11,gminsol21,gtmin1}=(<<solution1)[[{1,2,3}]];
		If[Max[sysDiff[0,T]]>tol,
			Print[""];difftmin;Print[notsatistmin];Print[SessionTime[]];
			{gminsol1,gminsol2,gtmin,gmindis,pmin,5epsilon}>>solution1,
			
			Print[""];difftmin;Print[satistmin];Print[SessionTime[]];
			maxErrorUptoT=Max[sysDiff[0,T]];tempT=IntegerPart[gtmin1];
			For[j=T+1,j≤IntegerPart[gtmin1],j++,If[Max[sysDiff[j-1,j]]>tol,tempT=j-1;Break[]];];
			Print[""];Print[step3];firstsol=0;precision=precision+4;increaseprec;
			{gminsol1,gminsol2,gtmin,gmindis,pmin,5epsilon}>>solution2;doneflag=0(* THIS DONEFLAG COMMAND SKIPS STEP 3 if set to 1 *);
			Goto[end]]],

	If[gtmin>((<<solution2)[[3]]-1)&& gmindis<(101/100)(<<solution2)[[4]],
	{gminsol11,gminsol21}=(<<solution2)[[{1,2}]];
	If[Max[sysDiff[0,T]]>tol,
		{gminsol1,gminsol2,gtmin,gmindis,pmin,epsilon}>>solution2;
		Print[""];diffprec;Print[notsatisprec];Print[SessionTime[]];precision=precision+4;increaseprec,

		{gminsol1,gminsol2,gtmin,gmindis,pmin,epsilon}>>solution1;
		Print[""];diffprec;Print[satisprec];Print[SessionTime[]];Print[""];Print[congratulations];Print[done];doneflag=1];
		Goto[end]];
		(* solution1 stores the final solution from step 3 and solution2 stores the final solution from step 2 unless there is a failure of precision in step 3. *);
	
];

rit++;Goto[begin];

Label[end];
{gminsol1,gminsol2,gtmin,gmindis,pmin,5epsilon}];

(* Set up for solving for the global nonlinear saddle path *)
solvesys[ini_,t1_,t2_]:=Block[{sol},

sol=NDSolve[{system,ini},Drop[vars,-1],{t,t1,t2},
		Method->{"EventLocator", "Event"->box,"EventAction":> Throw[Null, "StopIntegration"],"EventLocationMethod"->"StepBegin"},
			MaxSteps->Infinity,WorkingPrecision->digits,PrecisionGoal->Min[10,digits/2],AccuracyGoal->Infinity];
sol];

compdis[p_,ϵ_]:=Block[{solveit,theta,ntraj,pini,pexp,ini,sol1,sol2,dis,tdis,j,jdis,jtdis,dissol1,dissol2,jp,tab,mindis,mindis2,j1,j2,ii,tmin,jmin,minsol1,minsol2,tmin2,pmin1,jmin2,minsol12,minsol22,pmin2},

ntraj=ntrajini;firstit=1;

solveit:=(theta = 2π j/ntraj;pini=p+ϵ{Sin[theta], Cos[theta]};
If[Length[sysexpand]==0,pexp={},pexp=First[sysexpand][Apply[Sequence,Last[sysexpand]]];];

(* IF CLAUSE FOR OUT OF BOX SO THAT IT IS DONE ONLY FOR SOME INITIAL RINGS *)
If[rit≤5,
Table[vars[[i]][0]=.,{i,1,Length[vars]-1}];Evaluate[Table[vars[[i]][0],{i,1,Length[vars]-1}]]=Apply[Sequence,{pini,pexp,Drop[Drop[varso,2+Length[pexp]],-1]},1];
	If[(box/.t->0)== True,Print["Out of Box."];{dis,tdis}={(10^7),100};Goto[skip]];Table[vars[[i]][0]=.,{i,1,Length[vars]-1}];];

ini=Thread[Table[vars[[i]][0],{i,1,Length[vars]-1}]==Apply[Sequence,{pini,pexp,Drop[Drop[varso,2+Length[pexp]],-1]},1],List];
sol1=solvesys[ini,0,ts-(10^(-12))];If[sol1[[1,1,2,1,1,2]]<ts-(10^(-12)),{dis,tdis}=distance[sol1,sol2,sol1[[1,1,2,1,1,2]],0];];

If[sol1[[1,1,2,1,1,2]]== ts-(10^(-12)),
	ini = Thread[Table[vars[[i]][ts],{i,1,Length[vars]-1}]==Drop[Last[antishock][Apply[Sequence,Flatten[Table[vars[[i]][ts-(10^(-12))],{i,1,Length[vars]}]/.sol1]]],-1],List];
	sol2=solvesys[ini,ts,Tmax];{dis,tdis}=distance[sol1,sol2,sol2[[1,1,2,1,1,2]],If[firstit==1,0,Max[0,Min[tmin/2,tmin-5]]]];];Label[skip]);

If[ϵ==0,solveit;mindis2=dis;tmin2=tdis;pmin2=p;minsol12=sol1;minsol22=sol2;Goto[end2]];

For[j=1,j≤ ntraj  ,j++,
	If[trajtrack==0||(trajtrack==1 && (j==minpos[[1]]||j==minpos[[2]])),
		solveit;jdis[j]=dis;jtdis[j]=tdis;dissol1[j]=sol1;dissol2[j]=sol2;jp[j]=pini,jdis[j]=(10^7)];
];

firstit=0;
If[iflag==1 &&Length[minpos]==2,
trajtrack= 1,

tab=Table[jdis[Mod[ii-1,ntraj]+1]-jdis[Mod[ii,ntraj]+1],{ii,1,ntraj}];minpos=Flatten[Position[Table[tab[[Mod[ii-2,ntraj]+1]]>0 &&tab[[Mod[ii-1,ntraj]+1]]<0,{ii,1,ntraj}],True]];];

Print["Local minima from initial search over the ring."];
For[ii=1,
	ii≤Length[minpos],ii++,
	Print[ii,":" ,"   {mindis, tdis, jmin, ntraj} = ",{jdis[minpos[[ii]]],jtdis[minpos[[ii]]]//N,minpos[[ii]],ntraj}];
];

mindis=mindis2=(10^7);j1=1;j2= ntraj;

For[ii=1,ii≤Length[minpos],ii++,
	{mindis,tmin,jmin,ntraj,minsol1,minsol2}={jdis[minpos[[ii]]],jtdis[minpos[[ii]]],minpos[[ii]],ntrajini,dissol1[minpos[[ii]]],dissol2[minpos[[ii]]]};
	pmin1=jp[minpos[[ii]]];

	While[True,
	{j1,j2}=2{jmin-1,jmin+1};{ntraj,jmin}=2{ntraj,jmin};
	For[j=j1+1,j≤ j2  ,j+=2,
		solveit;
		If[ Max[Abs[Apply[(#1/#2-1)&,{pini,pmin1}]]]<(10^(-precision)),
			Print["You may be losing precision in calculations."] ;tmin2=(10^100);Goto[end2];];
		If[ Abs[mindis-dis]≤gmindis/100000,Goto[end1];];
		If[dis≤mindis,mindis=dis;tmin=tdis;jmin=j;minsol1=sol1;minsol2=sol2;pmin1=pini;Break[];];
	];
	];

	Label[end1];

	If[mindis<mindis2,
		mindis2=mindis;tmin2=tmin;jmin2=jmin;minsol12=minsol1;minsol22=minsol2;pmin2=pmin1];
	If[mindis2<gmindis ,Break[]];

];
Print["Ring ",rit,":    Minimum Distance = ",mindis2//N,"    Time of Closest Approach = ",tmin2//N,"    No. of Subdivisions of ϵ-ring > ",jmin2];

Label[end2];
{mindis2,tmin2,minsol12,minsol22,pmin2}];

sysDiff[t1_?NumberQ,t2_?NumberQ]:=Table[Max[Abs[Table[((vars[[i]][t]/.If[t<ts,gminsol1,gminsol2])-(vars[[i]][t]/.If[t<ts,gminsol11,gminsol21]))[[1]],{t,t1,t2,1/10}]]],{i,1,Length[vars]-1}];

durationCheck[___]:=Block[{j},
(* solution1 stores the final solution from step 3 and solution2 stores the final solution from step 2 unless there is a failure of precision in step 3. *)
{gminsol11,gminsol21}=(<<solution2)[[{1,2}]]; {gminsol1,gminsol2,gtmin}=(<<solution1)[[{1,2,3}]];

For[j=T+1,j≤IntegerPart[tempT],j++,
	If[Max[maxErrorUptoT,Max[sysDiff[j-1,j]]]>tol,
		Print["The solution is accurate within tolerance 'tol' = ",tol//N," upto ",j-1," years."];Break[]];];

If[j==IntegerPart[tempT]+1,Print["The solution is accurate within tolerance 'tol' = ",tol//N," upto ",j-1," years. Continue with next cell."];];
Print["The solution upto ",T," years is accurate to within ",Max[maxErrorUptoT,Max[sysDiff[0,T]]]//N,". Continue with next cell."];

];

End[];

(* Epilog *)
Protect[Evaluate[$Context<>"*"]];