epsilon=10^(-4);
exponent=-1;

delta=.1;
deltad=.1;
alpha=.67;
chi=1;
sigma=10;
eta=10;
phi=1;
beta=.99;
R=1/beta;
rho=.9;
thetac=1;
thetaf=1;

gamma=2;
A0=1;
kappa0=10^(-1);
M0=1;
    
sig2eps=.01;
sig2kappa=.01;

D=.66;



savedParameters=[epsilon exponent delta deltad alpha chi sigma eta phi beta R rho thetac thetaf gamma A0 kappa0 M0 sig2eps sig2kappa D];
save savedParameters
    
clear y
clear InitialValues
x0=ones(3,1);
options=optimset('maxIter',20000,'MaxFunEvals',20000, 'TolFun', 10^(-16));
InitialValues=fsolve(@System06_16_10,x0,options);
    
P=InitialValues(1);
lambda=InitialValues(2);
M=InitialValues(3);
xi=kappa0/(1-beta)/M;
    
C=(P*xi)^(-1/gamma);

N=((R-1+delta)*lambda/epsilon)^(1/exponent);             

Pd=(1-delta)*lambda*sigma*P/(sigma-1);
Pn=((P^(1-sigma)-D*Pd^(1-sigma))/(1-D))^(1/(1-sigma));
Ln=(C/A0*(Pn/P)^(-sigma))^(1/alpha);

w=(sigma-1)*P^(sigma-sigma/alpha)/sigma*alpha*A0^(1/alpha)*C^(1-1/alpha)*Pn^(1+sigma/alpha-sigma);
Ld=(w/P/A0/lambda/alpha)^(1/(alpha-1));
L=D*Ld+(1-D)*Ln;

V=P*C-w*L;
RE=(1-deltad)/deltad*V;
Div=V;
    

    
InitialValues=[P lambda M xi C N Pd Pn Ln w Ld L V RE Div];
save InitialValues;

dynare Attempt06_16_10.mod