%% steady state solver %%
%% Line 2 of Table 1 %%
clear all;
global U1;global U2;
U1=3.178; U2=2.281; U3=1.244; U4=1.064; U5=0.664;

x0 = 1;
options=optimset('Display','iter');   % Option to display output
[x,fval] = fsolve(@steadystate1,x0,options)  % Call solver

% prediction on U46 %
for i=1:47
    pu(i,1) = x.^i;
end

hat_U46 = round((U1*(sum(pu(3:5,1),1)))*1000); % report this with x 
hat_U712 = round((U1*(sum(pu(6:11,1),1)))*1000); % report this with x 
hat_U13 = round((U1*(sum(pu(12:47,1),1)))*1000); % report this with x 

table1 = zeros(9,10);
table1(1,6:10) = 1000*[U1,U2,U3,U4,U5];
table1(2,1:2) = [U1*1000,round(pu(1)*1000)/1000];
table1(2,6:10) = [U1*1000,U2*1000,hat_U46,hat_U712,hat_U13];
clearvars -except table1
%% Line 3 of Table 1 %%
global U1;global U2;global U3;global U4;
U1=3.178; U2=2.281; U3=1.244; U4=1.064; U5=0.664;

x0 = [2,1,1,1];
options=optimset('Display','iter');   % Option to display output
[x,fval] = fsolve(@steadystate3,x0,options)  % Call solver

wl = round(x(1)*1000);
wh = round(x(2)*1000);
pul = round(exp(-exp(x(3)))*1000)/1000;
puh = round(exp(-exp(x(4)))*1000)/1000;

for j=13:48
uhat(j,1) = x(1)*(pul^(j-1))+x(2)*(puh^(j-1));
end
hat_U13 = round((sum(uhat(13:48,1),1))*1000); % report this with x 

table1(3,1:4) = [wh,wl,puh,pul];
table1(3,6:10) = [1000*U1,1000*U2,1000*U3,1000*U4,hat_U13];
clearvars -except table1

%% Line 5 of Table 1 %%

global U1;global U2;global U3;global U4;global U5;
U1=3.178; U2=2.281; U3=1.244; U4=1.064; U5=0.664;

nn=1;
x0 = nn*[0.1;0.1;0.1;0.1;0.1];lb = [0.001;0.01;-5;-0.2;-0.5];ub =  [2;5;-0.01;3;1];

gs = GlobalSearch;
problem = createOptimProblem('fmincon','x0',x0,...
    'objective',@steadystate4,'Aineq',[],'bineq',[],'lb',lb,'ub',ub);
[x,fval,flag] = run(gs,problem);

wl = round(x(1)*1000);
wh = round(x(2)*1000);
pul = round(exp(-exp(x(3)))*1000)/1000;
puh = round(exp(-exp(x(4)))*1000)/1000;
d=round(x(5)*1000)/1000;

table1(5,1:5) = [wh,wl,puh,pul,d];
table1(5,6:10) = [1000*U1,1000*U2,1000*U3,1000*U4,1000*U5];
clearvars -except table1

%% Line 7 of Table 1 %%
global U1;global U2;global U3;global U4;
U1=3.339;U2=2.787;U3=2.131;U4=2.426;U5=1.902;
x0 = [2,1,1,1];
options=optimset('Display','iter');   % Option to display output
[x,fval] = fsolve(@steadystate3,x0,options)  % Call solver

% prediction on probs %
wl = round(x(1)*1000);
wh = round(x(2)*1000);
pul = round(exp(-exp(x(3)))*1000)/1000;
puh = round(exp(-exp(x(4)))*1000)/1000;

for j=13:48
uhat(j,1) = x(1)*(pul^(j-1))+x(2)*(puh^(j-1));
end

hat_U13 = round((sum(uhat(13:48,1),1))*1000); % report this with x 

table1(6,6:10) = 1000*[U1,U2,U3,U4,U5];
table1(7,1:4) = [wh,wl,puh,pul];
table1(7,6:10) = [1000*U1,1000*U2,1000*U3,1000*U4,hat_U13];
clearvars -except table1

%% Line 9 of Table 1 %%
global U1;global U2;global U3;global U4;global U5;
U1=3.339;U2=2.787;U3=2.131;U4=2.426;U5=1.902;

nn=1;
x0 = nn*[0.1;0.1;0.1;0.1;0.1];lb = [0.001;0.01;-5;-0.2;-0.5];ub =  [2;5;-0.01;3;1];

gs = GlobalSearch;
problem = createOptimProblem('fmincon','x0',x0,...
    'objective',@steadystate4,'Aineq',[],'bineq',[],'lb',lb,'ub',ub);
[x,fval,flag] = run(gs,problem);

wl = round(x(1)*1000);
wh = round(x(2)*1000);
pul = round(exp(-exp(x(3)))*1000)/1000;
puh = round(exp(-exp(x(4)))*1000)/1000;
d=round(x(5)*1000)/1000;

table1(9,1:5) = [wh,wl,puh,pul,d];
table1(9,6:10) = [1000*U1,1000*U2,1000*U3,1000*U4,1000*U5];
clearvars -except table1

%% Generates Line 4 and 8 of Table 1

Udata(1,:) = [3.178,2.281,1.244,1.064,0.664]; % Line 1: observed values
Udata(2,:) = [3.339,2.787,2.131,2.426,1.902]; % Line 6: observed values
 
for i=1:2
    global U1;global U2;global U3;global U4;global U5; global kout;global p;
    U1 = Udata(i,1);    U2 = Udata(i,2);    U3 = Udata(i,3);
    U4 = Udata(i,4);    U5 = Udata(i,5);
    
%% 5eqn/5unkown %%
x0 =[0.1;0.1;0.1;0.1];
lb =[-1;-1;-1;-1];
ub =[1;1;1;1];
kout = 0;
%gs = GlobalSearch;
%problem = createOptimProblem('fmincon','x0',x0,...
%    'objective',@steadystate4,'Aineq',[],'bineq',[],'lb',lb,'ub',ub);
%[x,fval,flag] = run(gs,problem);

[x,fval,flag] = fmincon(@steadystate_jh2,x0,[],[],[],[],lb,ub);
kout = 1;
steadystate_jh2(x);
kout = 0;

w = U1*1000;
d0=round(x(1)*1000)/1000;
d1=round(x(2)*1000)/1000;
d2=round(x(3)*1000)/1000;
d3=x(4);

nn = 4*i;
table1(nn,:) = [w,d0,d1,d2,d3,U1*1000,U2*1000,U3*1000,U4*1000,U5*1000];
P(:,i)=p;
end

