clear all
clc

%Parameters and options for later

b = 0.50;     
c = 0.50;     
rho = 0.50; 
rhox = 0.98;   
sigx = 0.014;    
rhoz = 0.70;    
sigz = 0.30;
beta = 0.95;    
alpha = 0.50;     
s = 0.50;     
xi = 0.97; % THIS CONTROLS RETURNS TO SCALE;     
delk = 0.10;     
delj = 0.30; 
gam0 = 10;
gam1 = -100;
order = 2; % ORDER FOR DYNARE SOLUTION;           
vcov = eye(2);

% finding steady state      
z = 0; % firm shocks turned off;                    
j = 0; % Inventory normalized to 1;
gam = gam0; % Price of risk const. in steady
p = log(delj*exp(j));    

% Solving for the remaining steady-state values;
InitGuess = [1,1];
tol = 12;
[SS_Values,fval] = fsolve(@(SS) steady(SS, delk, delj, c, b, alpha, rho, s, xi, beta),...
    InitGuess, optimset('TolX',10^(-tol),'TolFun',...
    10^(-tol),'MaxFunEvals',1000000,'MaxIter',50000,'algorithm','Levenberg-Marquardt'));
SS_Values = real(SS_Values);
k = SS_Values(1);
x = SS_Values(2);

% Finding the deterministic steady state of everything else;
i = log((((delk*exp(k)/c)^rho-(1-b)*((1-s)*exp(j))^rho)/b)^(1/rho));
CES = log(exp(k)-(1-delk)*exp(k));
CESi = log(c*(1/rho)*(b*exp(i)^rho+(1-b)*((1-s)*exp(j))^rho)^(1/rho-1)*b*rho*exp(i)^(rho-1) );
CESj = log(c*(1/rho)*(b*exp(i)^rho+(1-b)*((1-s)*exp(j))^rho)^(1/rho-1)*(1-b)*rho*((1-s)*exp(j))^(rho-1)*(1-s));
qk = log(1/exp(CESi));
xbar = x;
y = log(exp(x+z)*(exp(k)^alpha*(s*exp(j))^(1-alpha))^xi);
div = exp(y)-exp(i)-exp(p);
rf = log(1/beta); 

% Writing the mod file;
fid = fopen('stripped_model.mod','w+');
fprintf(fid,'var y m k i qk j p CES CESi CESj x z gam rf div; \n');
fprintf(fid,'varexo ex ez;\n');
fprintf(fid,'\n');
fprintf(fid, 'parameters delk delj alpha b c s rhoz sigz rhox sigx xbar gam0 gam1 beta rho xi;\n');
fprintf(fid,'delk = %3.12f;\n',delk);
fprintf(fid,'delj = %3.12f;\n',delj);
fprintf(fid,'c = %3.12f;\n',c);
fprintf(fid,'b = %3.12f;\n',b);
fprintf(fid, 'alpha = %3.12f;\n',alpha);
fprintf(fid,'s = %3.12f;\n',s);
fprintf(fid,'rhoz = %3.12f;\n',rhoz);
fprintf(fid,'sigz = %3.12f;\n',sigz);
fprintf(fid,'rhox = %3.12f;\n',rhox);
fprintf(fid,'sigx = %3.12f;\n',sigx);
fprintf(fid,'xbar = %3.12f;\n',xbar);
fprintf(fid,'gam0 = %3.12f;\n',gam0);
fprintf(fid,'gam1 = %3.12f;\n',gam1);
fprintf(fid,'beta = %3.12f;\n',beta);
fprintf(fid,'rho = %3.12f;\n',rho);
fprintf(fid, 'xi = %3.12f;\n',xi);
fprintf(fid, '\n');
fprintf(fid,'model;\n');
fprintf(fid,'exp(y) = exp(x+z)*(exp(k(-1))^alpha*(s*exp(j(-1)))^(1-alpha))^xi; // Production function \n');
fprintf(fid,'exp(j) = (1-delj)*exp(j(-1))+exp(p); // Accumulation of spare parts\n');
fprintf(fid,'exp(k) = (1-exp(delk))*exp(k(-1))+exp(CES); // Accumulation of physical capital \n');
fprintf(fid,'exp(CES) = c*(b*exp(i)^rho+(1-b)*((1-s)*exp(j))^rho)^(1/rho); // Adjustment cost (CES) \n');
fprintf(fid,'exp(CESi) = c*(1/rho)*(b*exp(i)^rho+(1-b)*((1-s)*exp(j))^rho)^(1/rho-1)*b*rho*exp(i)^(rho-1);\n');
fprintf(fid,'exp(CESj) = c*(1/rho)*(b*exp(i)^rho+(1-b)*((1-s)*exp(j))^rho)^(1/rho-1)*(1-b)*rho*((1-s)*exp(j))^(rho-1)*(1-s);\n');
fprintf(fid,'exp(qk) = 1/exp(CESi);\n');
fprintf(fid,'exp(qk) = exp(m(+1))*(xi*alpha*exp(y(+1))/(exp(k))+(1-exp(delk(+1)))*exp(qk(+1)));\n');
fprintf(fid,'1-exp(qk)*exp(CESj) = exp(m(+1))*(xi*(1-alpha)*exp(y(+1))/(exp(j))+1-delj);\n');
fprintf(fid,'exp(m) = exp(log(beta)+gam(-1)*(x(-1)-x));\n');
fprintf(fid,'gam = gam0 + gam1*(x-xbar);\n');
fprintf(fid,'z = rhoz*z(-1)+sigz*ez;\n');
fprintf(fid, 'x = (1-rhox)*xbar + rhox*x(-1)+sigx*ex;\n');
fprintf(fid,'1/exp(rf) = exp(m(+1));\n');
fprintf(fid,'div = exp(y)-exp(i)-exp(p);\n');
fprintf(fid,'end;\n');
fprintf(fid,'\n');

fprintf(fid,'initval;\n');
fprintf(fid,'y = %3.12f;\n',y);
fprintf(fid,'m = %3.12f;\n',log(beta));
fprintf(fid,'qk = %3.12f;\n',qk);
fprintf(fid,'k= %3.12f;\n',k);
fprintf(fid,'j= %3.12f;\n',k);
fprintf(fid,'p= %3.12f;\n',p);
fprintf(fid,'i= %3.12f;\n',i);
fprintf(fid,'x= %3.12f;\n',x);
fprintf(fid,'z= %3.12f;\n',z);
fprintf(fid,'gam= %3.12f;\n',gam);
fprintf(fid,'CES = %3.12f;\n',CES);
fprintf(fid,'CESi= %3.12f;\n',CESi);
fprintf(fid,'CESj= %3.12f;\n',CESj);
fprintf(fid,'rf= %3.12f;\n',rf);
fprintf(fid,'div= %3.12f;\n',div);
fprintf(fid,'end;\n');
fprintf(fid,'\n');

fprintf(fid, 'vcov =[ 1 0; 0 1;]; \n');
fprintf(fid, '\n');

fprintf(fid, 'order = %i;\n',order);
fclose(fid);
%%
% Trying to solve the model;
if exist('stripped_model.mat')
    delete stripped_model.mat 
end
!C:\dynare\4.5.0\dynare++\dynare++ --per 1 --sim 1 --irfs  ex ez stripped_model.mod