clear
clc

p.alpha1 = 0.5;
p.alpha2 = 1-p.alpha1;
p.alphan1 = 0.5; 
p.alphan2 = 1-p.alphan1;
p.rho = 0.9999; 
p.lambda = 1; 
p.eta = 1;
p.sigma = 1;  
p.beta = 1; 
p.epsilon = 6;

fun = @(x) root2d(x, p);
x0 = [1,0.5,0.5,1,1,1,1,0.5,0.5,1,0.5,0.5,1,1];
A = fsolve(fun,x0);
lb = zeros(size(x0)); %lower bound of zero
ub = Inf*ones(size(x0));
options = optimoptions('lsqnonlin','FunctionTolerance',1e-14,'OptimalityTolerance',1e-9);
[A,resnorm,residual,exitflag] = lsqnonlin(fun,x0,lb,ub,options)
B = [p.alpha1,p.alpha2,p.alphan1, p.alphan2, p.rho, p.lambda, p.eta, p.sigma, p.beta, p.epsilon];
A = [A,B];

save A A;

function F = root2d(x, p)
% hours and wages
F(1) = x(2)-p.alphan1*((x(5)/x(4))^p.lambda)*x(1);
F(2) = x(3)-p.alphan2*((x(6)/x(4))^p.lambda)*x(1);
F(3) = x(4)-(p.alphan1*(x(5)^(1+p.lambda))+p.alphan2*(x(6)^(1+p.lambda)))^(1/1+p.lambda);
% consumption
F(4) = x(8)-p.alpha1*x(7);
F(5) = x(9)-p.alpha2*x(7);
% labour supply
F(6) = x(4)-(x(1)^p.eta)*(x(7)^(1/p.sigma));
% production
F(7) = x(11)-x(13)*x(2);
F(8) = x(12)-x(14)*x(3);
% marginal costs
F(9) = (p.epsilon-1)/p.epsilon - x(5)/x(13);
F(10) = (p.epsilon-1)/p.epsilon - x(6)/x(14);
% equilibrium
F(11) = x(10)-x(11)-x(12);
F(12) = x(1)-x(2)-x(3);
F(13) = x(11)-x(8);
F(14) = x(12)-x(9);
end


