clear 
clc

%set parameters

r_annual = 6.5; %from K-R
rss = r_annual/400; % MAKE IT QUARTERLY
beta = 1/(1+rss);
delta = 0.021; %from G-W
gamma_a = 1.01; %assume a value
gamma_theta = 1.01; %assume a range of values
sigma = 0.15;
alpha = 1/3; 
rho = 0.95;
xi = 1; %similar to King-Rebello
M = 5; 


paramss = [sigma delta gamma_a gamma_theta alpha beta M];

z = ones(paramss(7),1);

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

z_steady = @(z)puttyclayss(paramss, z);

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

%compute the rest of the SS

for j = 1:M
    R(j) = beta^j * (gamma_a * gamma_theta^(alpha))^(j/(1-alpha)); 
    k_bar(j) = (1-delta)^(j-1) *R(j) * (1-normcdf(z(j)-sigma))*gamma_a^j;
    ylnum(j) = (1-delta)^(j-1) * (1-normcdf(z(j)-sigma))*(gamma_a*gamma_theta)^(-j*alpha/(1-alpha));
    ylden(j) = (1-delta)^(j-1) * (1-normcdf(z(j)));
    wnum(j) = (1-delta)^(j-1) *R(j) * (1-normcdf(z(j)-sigma))*gamma_a^j;
    wden(j) = (1-delta)^(j-1) *R(j) * (1-normcdf(z(j)))*(gamma_a*gamma_theta)^(j/(1-alpha));
    k_eff(j) = (1-delta)^(j-1) * (1-normcdf(z(j)))*(gamma_a*gamma_theta)^(-j/(1-alpha));
end 

kbar = (alpha *sum(k_bar))^(1/(1-alpha)); %this is capital embedded in a machine
xbar = kbar^alpha;
yl = sum(ylnum)/sum(ylden) * xbar;
L = 0.2; 
y = yl/L;
w = (1-alpha)* sum(wnum)/sum(wden) * xbar;
c = w*(1-L)/xi;
i = y - c;
Q = i/kbar;
ktilde = i* ((gamma_a * gamma_theta)^(1/(1-alpha)) -(1-delta))^(-1); %this is capital accumulation equation, unadjusted for endogenous depreciation
keff  = sum(k_eff)*kbar*Q; %this is the effective capital stock, where depreciation is an endogenous decision




