% Define parameters
chi = 1.5;  % substitutability of 
alpha = 0.33;  % capital share
beta = 0.99;  % discount rate
delta = 0.025;  % capital depreciation rate
psi = 0.4;  % positive externality of green-skilled labor
eta = 0.2;  % value of leisure

% Define gamma functions
gamma_n = @(j) 1;
gamma_g = @(j) 1.1845 * j.^2.12;

% Define the integrals with 'ArrayValued', true
integral_gamma_n = @(J) integral(@(j) (gamma_n(j)).^(chi-1), 0, J, 'ArrayValued', true);
integral_gamma_g = @(J) integral(@(j) (gamma_g(j)).^(chi-1), J, 1, 'ArrayValued', true);

% Define the function for aggregated labor services
L_function = @(J, L_n, L_g) ((integral_gamma_n(J)^(1/chi) * (L_n^(chi-1)/chi)) + ...
                             (integral_gamma_g(J)^(1/chi) * (L_g^(chi-1)/chi)))^(chi/(chi-1));

% Define the functions for W_n and W_g
W_n_function = @(L, L_n, J) (L / L_n)^(1/chi) * (integral_gamma_n(J)^(1/chi));
W_g_function = @(L, L_g, J) (L / L_g)^(1/chi) * (integral_gamma_g(J)^(1/chi));

% Define the equation for J
J_function = @(J, W_g, W_n) gamma_g(J) / gamma_n(J) - W_g / W_n;

% Define the function for Y
Y_function = @(K, L) K^alpha * L^(1-alpha);

% Define the functions for R and W
R_function = @(Y, K) alpha * (Y / K);
W_function = @(Y, L) (1 - alpha) * (Y / L);

% Household side conditions
Wg_household = @(C, L_g, L_n) C / (1 - L_g - L_n);
Wn_household = @(C, L_g, L_n) C / (1 - L_g - L_n);

% Euler equation
R_euler = 1 / beta - 1 + delta;

% Initial guesses
J0 = 0.5;  % Initial guess for J
L_n0 = 0.2;   % Initial guess for L_n
L_g0 = 0.2;   % Initial guess for L_g
Y0 = 1.2;     % Initial guess for Y
C0 = 0.9;     % Initial guess for C
K0 = 12;     % Initial guess for K

% Define the system of equations to solve
system_of_equations = @(x) [
    J_function(x(1), x(7), x(8));                    % Equation for J
    W_n_function(L_function(x(1), x(4), x(5)), x(4), x(1)) - x(8);  % Equation for W_n
    W_g_function(L_function(x(1), x(4), x(5)), x(5), x(1)) - x(7);  % Equation for W_g
    Y_function(x(3), L_function(x(1), x(4), x(5))) - x(2);          % Equation for Y
    Wg_household(x(6), x(5), x(4)) - x(7);           % Household condition for W_g
    Wn_household(x(6), x(5), x(4)) - x(8);           % Household condition for W_n
    x(2) - x(6) - delta * x(3);                      % Market clearing condition Y = C + delta * K
    R_function(x(2), x(3)) - R_euler;                % Euler equation for capital
    W_function(x(2), L_function(x(1), x(4), x(5))) - x(9)             % Wage definition
];

% Solve the system of equations
options = optimoptions('fsolve', 'Display', 'iter');
solution = fsolve(system_of_equations, [J0, Y0, K0, L_n0, L_g0, C0, W_g_function(L_function(J0, L_n0, L_g0), L_g0, J0), W_n_function(L_function(J0, L_n0, L_g0), L_n0, J0), W_function(Y0, L_function(J0, L_n0, L_g0))], options);

% Extract the solution
J_solution = solution(1);
Y_solution = solution(2);
K_solution = solution(3);
L_n_solution = solution(4);
L_g_solution = solution(5);
L_solution  = L_function(J_solution, L_n_solution, L_g_solution);
C_solution = solution(6);
W_g_solution = solution(7);
W_n_solution = solution(8);
W_solution = solution(9);

% Display the results
fprintf('J: %f\n', J_solution);
fprintf('Y: %f\n', Y_solution);
fprintf('K: %f\n', K_solution);
fprintf('L_n: %f\n', L_n_solution);
fprintf('L_g: %f\n', L_g_solution);
fprintf('C: %f\n', C_solution);
fprintf('W_g: %f\n', W_g_solution);
fprintf('W_n: %f\n', W_n_solution);
fprintf('L: %f\n', L_solution);
fprintf('W: %f\n', W_solution);
fprintf('R: %f\n', R_euler);