clear 
clc 

sigma = 0.5;
alpha = 1/3; 
r = 0.065;
beta = 1/(1+r);
psi = 1;
delta = 0.025;
rho = 0.95;
g = 0;
theta_ss = 1;  %(1)
M = 3;


params = [sigma alpha r beta psi delta rho g M];

z_ss = zeros(params(9),1);

options=optimset('display','iter');%,'TolFun',1e-4,'TolX', 1e-4);
options=optimset(options,'MaxFunEvals',10000,'TolX',1e-8);

z_steady = @(z_ss)steadygrow(params, z_ss);

[z_ss, fval1, exitflag1] = fsolve(z_steady, z_ss, options); %(2)
exitflag1

zssCell = num2cell(z_ss);
for j = 1:M
eval(sprintf('zss_%d = z_ss(j)', j))
end 

k1 = zeros(params(9),1);
y1 = zeros(params(9),1);
kl1 = zeros(params(9),1);
l1 = zeros(params(9),1);
for i = 1:params(9)
    k1(i) = (1-normcdf(z_ss(i)-params(1)))*(params(4)*(1-params(6)))^(i-1)*(1+params(8))^(-i);
    y1(i) = (1-normcdf(z_ss(i)-params(1)))*(1-params(6))^(i-1)*(1+params(8))^(-i);
    kl1(i) = (1-normcdf(z_ss(i)))*(params(4)*(1-params(6)))^(i-1);
    l1(i) = (1-delta)^(i-1)*(1-normcdf(z_ss(i)));
end
kss = (alpha* theta_ss *sum(k1))^(1/(1-alpha));
xss = theta_ss * kss^(alpha);
ratiok = sum(k1)/sum(kl1);
wss = (1-alpha) * ratiok*xss;
ratioy = sum(y1)/sum(l1);
y_lss = ratioy*xss;
q_lss = sum(l1)^(-1);
c_lss = y_lss - kss*q_lss;
lss = wss * (psi*c_lss)^(-1) * (1 + (wss * (psi*c_lss)^(-1)))^(-1);
yss = y_lss * lss;
css = c_lss * lss;
qss = q_lss * lss;

checkzz = zeros(params(9),1);
checkz = zeros(params(9),1);
for i = 1:params(9)
checkzz(i) = z_ss(i) - 1/sigma*((log(wss) - log (xss) +1/2 *sigma^2) + i * log(1+g));
checkz(i) = z_ss(i) - 1/sigma * (log(1-alpha) + log(sum(k1)) - log(sum(kl1)) + 1/2 * sigma^2 +i * log(1+g));
end 
checkzz
checkz

check = zeros(8,1);

check(1) = yss -sum(y1)*xss*qss;
check(2) = kss -  alpha * xss * sum(k1);
check(3) = wss - psi * css/(1-lss);
check(4) = css - yss + kss * qss;
check(5) = kss - (xss * sum(k1) - wss * sum(kl1));
check(6) = lss - qss * sum(l1);
check(7) = xss - theta_ss * kss^alpha;
check(8) = theta_ss -1;

vars = [theta_ss kss xss wss yss qss css lss]'

check

save parameterfile.mat

