using TwoStageOptimalControl

U(Con, Stat, t::Float64, Para::Dict) = Stat[2] * ((Para["A"] - Stat[1])Stat[1] - Para["b"]Con[1] - Para["c"]/2Con[1]^2) Q(Con,Stat, t::Float64,s::Float64, Para::Dict) = Stat[2] * ((Para["A"] - Stat[1])Stat[1] - Para["b"]Con[1] - Para["c"]/2Con[1]^2)

f1(Con,Stat, t::Float64, Para::Dict) = [Con[1] - Para["delta"]Stat[1], -Para["eta"]Stat[2]] f2(Con,Stat, t::Float64, s::Float64, Para::Dict) = [Con[1] - Para["delta"]*Stat[1], 0]

g(Con,Stat,t::Float64, Para::Dict) = [Stat[1], Para["eta"]*Stat[2]]

S1(Stat, Para::Dict) = 0 S2(Stat_dist, Para::Dict) = 0

MyPara = Dict() MyPara["A"] = 60 MyPara["b"] = 100 MyPara["c"] = 10 MyPara["delta"] = 0.1 MyPara["rho"] = 0.1 MyPara["eta"] = 0.05

MyPara["T"] = 10 MyPara["hstep"] = 0.5 MyPara["nCon"] = 1 MyPara["nCondist"] = 1 MyPara["nStat"] = 2 MyPara["nStatdist"] = 2 MyPara["InitStat"] = [0.1,1.0] MyPara["ConMin"] = [0.0] MyPara["Con_distMin"] = [0.0]

MyPara["OptiType"] = "Newton-Raphson" MyPara["ProbIndex"] = 2

MyPara["InitLineStep"] = 1e-5 MyPara["UpperLineStep"] = 1e-2 MyPara["hLowBound"] = 0.1 MyPara["PlotResultsIntermediateFrequency"] = 150

MyPara["LoadInits"] = true

Results = Dict() Results["Con"] = 6.0ones(1,10,1) Results["Con_dist"] = 6.0ones(10,10,1) Results = TwoStageOptimisation(Results = Results,UserParameters = MyPara, ObjectiveIntegrand = U, AggregationFunction = Q, StateDynamic1 = f1, StateDynamic2 = f2, Shock = g, SalvageFunction1=S1, SalvageFunction2=S2)

SaveResults(Results,"test/CapitalAccumulationBenchmark")